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Abstract 

We  analyze  bathymetric  and  gravity  anomalies  at  five  plume-ridge  systems  to  constrain 
cmstal  and  mantle  density  structure  at  these  prominent  oceanic  features.  Numerical  models 
are  then  used  to  explore  the  physical  mechanisms  controlling  plume-ridge  interaction  and  to 
place  theoretical  constraints  on  the  temperature  anomalies,  dimensions,  and  fluxes  of  the 
Icelandic  and  Galapagos  plumes. 

In  Chapter  1  we  analyze  bathymetric  and  gravity  anomalies  along  the  hotspot- 
influenced  Galapagos  Spreading  Center.  We  find  that  the  Galapagos  plume  generates 
along-axis  bathymetric  and  mantle-Bouguer  gravity  anomalies  (MBA)  that  extend  >500  km 
east  and  west  of  the  Galapagos  Islands.  The  along-axis  MBA  becomes  increasingly 
negative  towards  the  plume  center,  reaching  a  minimum  of  — 90  mGal  near  91°W,  and  axial 
topography  shallows  by  -1.1  km  toward  the  plume.  These  variations  in  MBA  and 
bathymetry  are  attributed  to  the  combined  effects  of  cmstal  thickening  and  anomalously  low 
mantle  densities,  both  of  which  are  due  to  a  mantle  temperature  anomaly  imposed  beneath 
the  ridge  by  the  Galapagos  plume.  Passive  mantle  flow  models  predict  a  temperature 
anomaly  of  50±25°C  is  sufficient  to  produce  the  2-4  km  excess  cmst  required  to  explain  the 
along-axis  anomalies.  70-75%  of  the  along-axis  bathymetric  and  MBA  variations  are 
estimated  to  arise  from  the  cmst  with  the  remaining  25-30%  generated  by  the  anomalously 
hot,  thus  low-density  mantle.  Along  Cocos-plate  isochrons,  bathymetric  and  MBA 
variations  increase  with  increasing  isochron  age,  suggesting  the  subaxial  mantle 
temperature  anomaly  was  greater  in  the  past  when  the  plume  was  closer  to  the  ridge  axis. 

In  addition  to  the  Galapagos  plume-ridge  system,  in  Chapter  2  we  examine  along- 
isochron  bathymetric  and  MBA  variations  at  four  other  plume-ridge  systems  associated 
with  the  Iceland,  Azores,  Easter  and  Tristan  hotspots.  We  show  that  residual  bathymetry 
(up  to  4.7  km)  and  mantle-Bouguer  gravity  anomalies  (up  to  -340  mGal)  are  greatest  at  on- 
axis  plumes  and  decreases  with  increasing  ridge-hotspot  separation  distance,  until 
becoming  insignificant  at  a  plume-ridge  separation  of  -500  km.  Along-isochron  widths  of 
bathymetric  anomalies  (up  to  2700  km)  decrease  with  increasing  paleo-spreading  rate, 
reflecting  the  extent  to  which  plume  material  flows  along-axis  before  being  swept  away  by 
the  spreading  lithosphere.  Scaling  arguments  suggest  an  average  ridgeward  plume  flux  of 
-2.2x1 0^  km/my.  Assuming  that  the  amplitudes  of  the  MBA  and  bathymetric  anomalies 
reflect  crustal  thickness  and  mantle  density  variations,  passive  mantle  flow  models  predict 
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maximum  subaxial  mantle  temperature  anomalies  to  be  150-225°C  for  ridge-center  plumes, 
which  decrease  as  the  ridges  migrate  away  from  the  plumes. 

The  dynamics  of  mantle  flow  and  melting  at  ridge-centered  plumes  are  investigated  in 
Chapters  3  using  three-dimensional,  variable-viscosity,  numerical  models.  Three 
buoyancy  sources  are  examined:  temperature,  melt  depletion,  and  melt  retention.  The 
width  W  to  which  a  plume  spreads  along  a  ridge  axis  depends  on  plume  volume  flux  Q, 
full  spreading  rate  U,  buoyancy  number  B  =  (QApg)/(4SrioU'^),  and  ambient/plume 
viscosity  contrast  /according  to  lT=2.37((2/f/) ^^^(57)0.04  Thermal  buoyancy  is  first 
order  in  controlling  along-axis  plume  spreading  while  latent  heat  loss  due  to  melting,  and 
depletion  and  retention  buoyancy  forces  contribute  second  order  effects.  Two  end-member 
models  of  the  Iceland-Mid-Atlantic  Ridge  (MAR)  system  are  examined.  The  first  end- 
member  model  has  a  broad  plume  source  of  radius  300  km,  temperature  anomaly  of  75°C, 
and  volume  flux  of  1.2x10'  km^/my.  The  second  model  has  a  narrower  plume  source  of 
radius  60  km,  temperature  anomaly  of  170°C,  and  flux  of  2.1x10^  km^/my.  The  first 
model  predicts  successfully  the  observed  crustal  thickness,  topographic,  and  MBA 
variations  along  the  MAR,  but  the  second  model  requires  substantial  along-axis  melt 
transport  in  order  to  explain  the  observed  along-axis  variations  in  crustal  thickness, 
bathymetry,  and  gravity.  We  favor  this  second  model  because  it  predicts  a  mantle  P-wave 
velocity  reduction  in  the  plume  of  ~2%  as  consistent  with  recent  seismic  observations 
beneath  Iceland. 

Finally  in  Chapter  4  we  use  three-dimensional  numerical  models  to  investigate  the 
interaction  of  plumes  and  migrating  midocean  ridges.  Scaling  laws  of  axial  plume 
spreading  width  W  are  derived  first  for  stationary  ridges  and  off-axis  plumes,  which  yield 
results  consistent  with  those  obtained  from  independent  studies  of  Ribe  [1996].  W  and  the 
maximum  plume-ridge  interaction  distance X/nox  again  scale  with  {QIU)^^  as  in  the  case  of 
ridge-centered  plumes  and  increase  with  7  and  buoyancy  number.  In  the  case  of  a 
migrating  ridge,  x^ax  is  reduced  when  a  ridge  migrates  toward  the  plume  due  to  excess 
drag  of  the  faster-moving  leading  plate,  and  enhanced  when  a  ridge  migrates  away  from  the 
plume  due  to  reduced  drag  of  the  slower-moving  trailing  plate.  Thermal  erosion  of  the 
lithospheric  boundary  layer  by  the  previously  ridge-centered  plume  further  enhances  W  and 
Xmax  but  to  a  degree  that  is  secondary  to  the  differential  migration  rates  of  the  two  plates. 
Model  predictions  are  compared  with  observed  along-isochron  bathymetric  and  MBA 
variations  at  the  Galapagos  plume-ridge  system.  The  anomaly  amplitudes  and  widths,  as 
well  as  the  increase  in  anomaly  amplitude  with  age  are  predicted  with  a  plume  source 
temperature  anomaly  of  80-120°C,  radius  of  80-100  km,  and  volume  flux  of  4.5x10^ 
km^'m.y.  Our  numerical  models  also  predict  crustal  production  rates  of  the  Galapagos 
Islands  consistent  with  those  estimated  independently  using  the  observed  island 
topography.  Predictions  of  the  geochemical  signature  of  the  plume  along  the  present-day 
ridge  suggest  that  mixing  between  the  plume  and  ambient  mantle  sources  is  unlikely  to 
occur  in  the  asthenosphere  or  shallow  crust,  but  most  likely  deeper  in  the  mantle  possibly 
by  entrainment  of  ambient  mantle  as  the  plume  ascends  through  the  depleted  portion  of  the 
mantle  from  its  deep  source  reservoir. 

Thesis  Supervisor:  Dr.  Jian  Lin 

Department  of  Geology  and  Geophysics 
Wood  Hole  Oceanographic  Institution 
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Introduction 


Hotspots  and  midocean  ridges  are  the  sources  of  the  ocean's  igneous  crust  and  are  thus 
the  primary  mechanisms  by  which  heat  and  mass  are  transported  from  the  mantle  to  the 
Earth's  surface.  The  present-day  global  oceanic  crustal  production  rate  of  2  x  10^  km^/my, 
of  which  -10%  is  due  to  hotspot  volcanism  [Larson  ,  1991]  is  sufficient  to  resurface  the 
planet  with  a  7  km-thick  crust  every  175  m.y.  years.  Moreover,  crustal  production  rates 
may  have  been  greater  by  a  factor  of  2  in  the  geologically  recent  past  [Larson,  1991]. 
Thus,  studies  of  igneous  and  mantle  dynamic  processes  at  hotspots  and  midocean  ridges 
are  crucial  to  our  understanding  of  Earth  structure  at  present-day  and  in  the  past. 

Over  the  past  three  decades  much  has  been  learned  about  the  dynamics  of  mantle  flow 
and  melt  generation  at  hotspots  and  midocean  ridges.  Since  Hess's  [1962]  hypothesis  that 
midocean  ridges  are  the  ascending  limbs  of  mantle  convection  cells,  a  number  of 
observational  and  theoretical  studies  have  shaped  our  present  conceptions  of  midocean 
ridge  dynamics.  For  example,  Hess's  [1962]  convection  hypothesis  was  examined 
quantitatively  by  means  of  a  boundary  layer  treatment  of  cell  convection  in  two-dimensions 
(2-D)  [Oxburgh  and  Turcotte,  1967;  Turcotte  and  Oxburgh,  1967].  The  study  by  Oxburgh 
and  Turcotte  [1967]  was  among  the  first  to  establish  the  concept  of  a  lithospheric  thermal 
boundary  layer,  to  explain  the  decrease  in  seafloor  heat  flow  with  age,  and  to  discuss 
decompression  melting  processes  at  midocean  ridges.  A  parallel  study  by  McKenzie 
[1967]  was  among  the  first  to  explain  seafloor  heat  flow  variations  by  a  conductively 
cooling  plate  model  overlying  an  asthenosphere  of  uniform  temperature. 

Furthermore,  the  finding  that  normal  oceanic  crust  was  -6  km  in  thickness,  globally 
[Raitt,  1963] ,  provides  a  powerful  constraint  on  mantle  flow  and  thermal  structure  beneath 
midocean  ridges.  Reid  and  Jackson  [1981]  demonstrated  that  simple  2-D  corner  flow 
models  could  produce  the  mantle  temperatures  and  upwelling  rates  necessary  to  generate  a 
6-km  thick  crust  at  intermediate  and  fast  spreading  rates.  Thermal  and  compositional 
buoyancy,  however,  seemed  to  be  required  to  generate  6  km  of  crust  at  slow  spreading 
rates  [e.g.  Buck  and  Su,  1989;  Rabinowicz,  1987;  Scott  and  Stevenson,  1989;  Sotin  and 
Parmentier,  1989].  Further  work  such  as  that  of  Bottinga  andAllegre  [1978],  Klein  and 
Langmuir  [1987],  McKenzie  [1984],  md  McKenzie  and  Bickle,  [1988]  were  groundwork 
studies  of  the  thermal  dynamics  of  mantle  melting  at  ridges  and  on  the  composition  of 
ocean  ridge  basalts. 
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During  the  time  that  seafloor  spreading  was  being  recognized  as  the  surface  expression 
of  a  convecting  Earth,  hotspots  were  also  being  attributed  to  mantle  processes,  separate 
from,  but  consistent  with,  the  plate  tectonic  paradigm.  Wilson  [1963]  introduced  the 
concept  that  the  age  progression  along  the  Hawaiian  Island  chain  reflects  migration  of  the 
lithospheric  plate  over  a  magma  source  which  is  fixed  in  the  mantle.  Such  findings  led  to 
Morgan's  hypothesis  that  hotspots  are  the  result  of  mantle  convection  plumes  which  ascend 
from  the  deep  mantle  to  the  base  of  the  lithosphere  [Morgan,  1971;  Morgan,  1972]. 
Follow-up  theoretical  studies  by  Parmentier  et  al.  [1975]  and  laboratory  experiments  of 
Whitehead  and  Luther  [1915]  demonstrated  that  Rayleigh-Taylor  instabilities  from  a  basally 
heated  mantle  could  generate  localized  upwelling  in  the  form  of  plumes.  Along  with  the 
plume  hypothesis  came  studies  by  Crough  [1978,  1983],  who  attributed  the  broad  (1000- 
2000  km)  hotspot  swells  to  anomalously  hot,  low-density  lithosphere,  and  work  by  Detrick 
and  Crough  [1978],  which  introduced  the  concept  of  plume-induced  lithospheric  thinning. 
Later  work  by  Olson  and  colleagues  showed  that  hotspot  swells  could  be  explained  by  the 
dynamic  uplift  of  a  plume  as  it  spreads  gravitationally  beneath  the  lithosphere  [Olson  1990; 
Olson  and  Nam,  1986;  Olson  et  al.,  1988].  Finally,  Watson  and  McKenzie  [1991] 
combined  the  physics  of  a  buoyantly  upwelling  plume  with  melting  models  of  McKenzie 
and  Bickle  [1988]  to  examine  melting  processes  beneath  the  Hawaiian  hotspot. 

A  landmark  discovery  by  Schilling  and  co-workers  demonstrated  that  igneous  products 
at  hotspots  such  as  Hawaii,  Iceland,  Galapagos,  and  the  Azores  have  rare-earth  element 
compositions  distinct  from  typical  midocean  ridge  basalts  (MORE)  [Schilling,  1971,  1973, 
1975,  Schilling  et  al,  1976;  Schilling  and  Winchester,  1967].  Moreover,  their  findings  of 
hotspot-type  chemical  signatures  in  basalts  at  hotspot-like  swells  along  midocean  ridges, 
such  as  Iceland  on  the  Mid-Atlantic  Ridge,  led  to  the  concept  that  rising  mantle  plumes 
interact  with  and  feed  oceanic  spreading  centers  [Hartetal,  1973;  Schilling,  1971,  1973, 
1975;  Schilling  and  Winchester,  1967;  Sun  et  al,  1975].  Independent  studies  of  plate 
kinematics  led  to  Morgan  [1978]'s  idea  of  a  "second  type  of  hotspot  island"  which  also 
suggested  that  plumes  spread  horizontally  to  nearby  oceanic  spreading  centers;  while  Vogt 
[1971,  1972,  1976]  showed  evidence  that  plumes  inject  material  also  along  the  axes  of 
midocean  ridges.  These  original  studies  stimulated  numerous  geophysical  and  geochemical 
surveys  of  plume-ridge  systems  leading  to  papers  by  Schilling  and  co-workers  which  have 
shaped  concepts  today  of  how  mantle  plumes  may  interact  with  midocean  ridges  [e.g. 
Schilling,  1985,  1991;  Schilling  et  al,  1985] . 
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With  the  conceptual  frameworks  of  midocean  ridges,  plumes,  and  the  interaction  of 
plumes  and  ridges  established,  the  purpose  of  this  thesis  is  to  examine  quantitatively  the 
mantle  and  crustal  structure  of  plume-ridge  systems  and  the  causal  mantle  dynamic  and 
igneous  processes.  The  first  two  chapters  focus  on  using  bathymetric  and  gravity 
observations  to  infer  crustal  and  mantle  density  structure  at  the  Galapagos  (Chapter  1), 
Iceland,  Azores,  Tristan,  and  Easter  (Chapter  2)  plume-ridge  systems.  The  last  two 
chapters  focus  on  the  dynamics  of  mantle  flow  and  melting  at  plume-ridge  systems,  which 
are  investigated  with  numerical  models  as  constrained  by  the  geophysical  observations. 

In  Chapter  1,  we  investigate  the  crustal  thickness  and  mantle  temperature  variations 
along  the  Galapagos  Spreading  Center  imposed  by  the  Galapagos  plume.  The  mantle- 
Bouguer  gravity  anomaly  (MBA) — which  is  the  free-airy  gravity  anomaly  corrected  for  the 
attraction  of  seafloor  topography  and  the  crust-mantle  interface  assuming  a  reference  crust 
of  uniform  density  and  thickness — has  been  particularly  useful  in  understanding  subsurface 
density  structure  at  midocean  ridges.  For  example,  bulls-eye  shaped  MBA  lows  centered 
on  individual  ridge  segments  as  documented  by  Kuo  and  Forsyth  [1988]  and  Lin  et  al. 
[1990]  are  strong  evidence  that  crustal  accretion  at  slow-spreading  ridges  varies 
significantly  along-axis  and  that  this  accretion  may  occur  due  to  convective  upwelling  as 
hypothesized  by  Whitehead  et  al.  [1984].  In  Chapter  1  we  produce  maps  of  MBA  anomaly 
and  bathymetry,  both  of  which  reflect  variations  in  crustal  thickness  and  mantle  density  at 
the  Galapagos  ridge  due  to  the  excess  temperature  imposed  by  the  Galapagos  hotspot. 

Temperature  anomalies  and  the  structure  of  mantle  plumes  at  intraplate  hotspots  are 
reflected  directly  in  the  shape  and  amplitude  of  hotspot  swells  [e.g.  McNutt,  1987;  Sleep 
1987,  1990].  At  near-ridge  hotspots  such  as  the  Galapagos,  however,  the  mantle 
temperature  anomaly  at  the  ridge-axis  is  likely  to  enhance  crustal  production;  consequently, 
investigations  of  mantle  temperature  anomalies  at  near-ridge  hotspots  requires  the 
consideration  of  melting.  We  therefore  use  3-D  passive  flow  models — analogous  to  the 
comer  flow  description  of  Reid  and  Jackson  [1981] — of  the  Galapagos  Spreading  Center 
to  predict  thickness  variations  along  the  ridge  axis  due  to  a  range  of  axial  temperature 
anomalies.  The  combined  contributions  of  crustal  thickness  and  mantle  density  variations 
to  bathymetry  and  MBA  are  then  compared  with  observations  to  constrain  mantle 
temperature.  Similar  analyses  are  done  for  anomalies  along  Cocos-Plate  isochrons  to  infer 
cmstal  thickness  variations  and  associated  mantle  temperature  anomalies  in  the  past  when 
the  Galapagos  Spreading  Center  was  closer  to  the  Galapagos  plume. 


13 


Analyses  of  bathymetric  and  MBA  variations  along  isochrons  is  a  unique  method  of 
investigating  the  evolution  of  individual  plume-ridge  systems.  In  Chapter  2  we  investigate 
along-axis  and  along-isochron  anomalies  at  five  prominent  systems  associated  with  the 
Galapagos,  Azores,  Iceland,  Tristan,  and  Easter  hotspots.  In  addition  we  use  independent 
constraints  on  past  plate  motions  to  estimate  plume-ridge  separation  distances  and 
spreading  rates  at  times  corresponding  to  the  isochron  ages.  We  investigate  relationships 
between  bathymetric  and  MBA  amplitudes  and  paleo-plume-ridge  distance,  as  well  as 
between  widths  of  along-isochron  bathymetric  anomalies  and  paleo-spreading  rate.  Scaling 
laws  are  then  derived  for  the  dependence  of  anomaly  amplitudes  and  mantle  temperature 
anomalies  to  examine  how  axial  temperature  anomalies  of  the  five  systems  may  have 
changed  with  plume-ridge  distance. 

While  the  passive  flow  models  used  in  Chapters  1  and  2  are  reasonable  approximations 
of  the  flow  beneath  oceanic  spreading  centers,  they  are  poor  representations  of  the  flow 
structure  at  buoyantly  upwelling  plumes.  To  investigate  the  dynamics  of  mantle  flow  and 
melting  at  plume-ridge  systems  it  is  therefore  necessary  to  incorporate  both  the  flow 
beneath  a  spreading  center  system  as  first  investigated  by  Turcotte  and  Oxburgh,  as  well  as 
the  pertinent  physics  of  plume  convection  as  originally  identified  by  Parmentier  et  al., 
[1975]  and  Whitehead  and  Luther,  [1975] .  Thus,  in  Chapter  3  we  used  numerical  models 
to  simulate  the  3-D  interaction  between  ridge-centered  buoyant  plumes  and  oceanic 
spreading  centers.  We  consider  fully  pressure-  and  temperature-dependent  rheology  and 
investigate  buoyancy  due  to  thermal  expansion,  melt  depletion,  and  melt  retention.  First, 
scaling  laws  tire  derived  for  the  dependence  of  along-axis  plume  width  on  plume  flux,  ridge 
spreading  rate,  and  ambient/plume  viscosity  contrast  in  the  absence  of  melting.  We  then 
investigate  the  melting  effects  of  latent  heat  loss,  and  depletion  and  retention  buoyancy  on 
flow  structure  and  on  the  scaling  laws.  Finally,  we  apply  our  models  to  the  Iceland-Mid- 
Atlantic  Ridge  system.  Model  predictions  and  observations  of  along-axis  crustal  thickness, 
bathymetry,  MBA,  and  geochemical  variations  are  compared  for  two  plume  source  radii 
and  temperature  anomalies  which  represent  end-member  properties  of  the  Icelandic  plume 
source. 

The  purpose  of  the  last  chapter  is  to  investigate  the  fluid  dynamics  of  plume-migrating 
ridge  interaction.  An  important  aspect  is  to  test  quantitatively  the  "mantle  plume 
source/migrating  ridge  sink"  hypothesis  originally  proposed  by  Schilling  and  co-workers 
as  based  on  their  geochemical  findings  as  well  as  the  work  of  Morgan  and  Vogt.  Scaling 
laws  are  first  derived  for  off-axis  plumes  in  steady  state  with  stationary  midocean  ridges 
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and  are  compared  with  independent  but  parallel  studies  of  Ribe  [1996].  For  migrating 
ridge  cases,  we  then  investigate  how  along-axis  plume  width  and  maximum  plume-ridge 
interaction  distance  scale  with  ridge  migration  rate,  spreading  rate,  plume  flux,  and 
ambient/plume  viscosity.  The  model  is  then  tested  by  comparing  model  predictions  with 
bathymetric  and  MBA  observations  of  the  Galapagos  plume-migrating  ridge  system. 
Anomalies  are  compared  at  the  present-day  ridge  axis  as  well  as  at  the  Cocos  Plate 
isochrons  examined  in  Chapter  1.  We  also  compare  predictions  and  observations  of  the 
geochemical  signature  of  the  Galapagos  plume  along  the  Galapagos  Spreading  Center  in 
order  to  investigate  the  processes  of  mixing  between  the  plume  and  ambient  mantle 
sources. 

Finally  I  include  in  the  Appendix,  laboratory  tank  experiments  done  with  C.  Kincaid 
and  C.  Gable  on  off-axis  plume-ridge  interaction.  A  plume  and  ridge  upper  mantle  system 
is  simulated  with  a  tank  of  concentrated  sucrose  solution  in  order  to  investigate  the  primary 
mechanisms  that  allow  an  off-axis  plume  to  overcome  the  viscous  drag  of  a  spreading  plate 
to  feed  the  nearby  ridge. 
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Abstract  To  better  understand  the  effects  of  hot  spots  on  mid-ocean  ridge  thermal  structure,  we 
investigate  the  subsurface  density  structure  of  the  Galdpagos  spreading  center  and  nearby 
lithosphere.  Using  shipboard  gravity  and  bathymetry  data,  we  obtain  maps  of  mantle  Bouguer 
anomalies  (MBA)  by  removing  from  the  free-air  gravity  the  attractions  of  seafloor  topography 
and  a  6-km-thick  model  crust.  Comparison  of  observed  and  theoretical  MBA  profiles  along 
isochrons  for  ages  0.0-7.7  Ma  suggests  that  seafloor  topography  is  isostatically  compensated  by 
mass  anomalies  primarily  in  the  upper  100  km  of  the  mantle.  This  result  is  consistent  with  the 
notion  that  seafloor  topography  along  the  Galapagos  spreading  center  is  supported  by  lateral 
changes  of  crustal  thickness  and  upper  mantle  density,  both  of  which  are  controlled  by 
temperatures  in  the  upper  mantle  where  decompression  melting  occurs.  Along  the  ridge  axis,  the 
MBA  decreases  from  the  east  and  west  toward  the  Galapagos  hot  spot  by  -90  mGal,  reaching  a 
minimum  nearest  the  hot  spot  at  91°W.  Seafloor  topography  mirrors  the  MBA  along  axis, 
increasing  by  -1. 1  km  toward  the  hot  spot.  These  variations  in  MBA  and  bathymetry  can  be 
explained  by  crustal  thickening  and  mantle  density  variations  resulting  from  a  gradual  axial 
temperature  increase  of  50±25°C  toward  the  hot  spot.  The  predicted  crustal  thickening  of  2-4 
km  nearest  the  hot  spot  accounts  for  70-75%  of  the  along-axis  MBA  and  bathymetry  anomalies; 
mantle  density  variations  account  for  the  rest  of  the  anomalies.  From  the  crustal  isochron  of  age 
7.7  Ma  to  the  present-day  axis,  the  along-isochron  amplitudes  of  MBA  decrease  from  -150  to 
-90  mGal.  The  corresponding  along-isochron  bathymetry  anomalies  decrease  from  -1.7  to  -1.1 
km.  These  observations  along  the  paleoaxes  of  the  Galapagos  spreading  center  indicate  that  the 
axial  temperature  anomaly  was  70%  hotter  in  the  past  (86±25°C)  and  has  steadily  decreased  to 
50±25°C  as  the  ridge  axis  migrated  away  from  the  Galapagos  hot  spot.  These  along-isochron 
temperature  anomalies,  however,  have  remained  well  below  that  estimated  for  the  hot  spot  itself 
(200°C),  indicating  that  the  lateral  temperature  gradient  between  the  hot  spot  and  the  ridge  axis 
has  remained  10-20  times  greater  than  that  along  the  ridge  axis  over  the  past  7.7  m.y. 


Introduction 

Three-dimensional  gravity  studies  of  mid-ocean  spreading 
centers  have  proven  crucial  to  understanding  the  processes 
controlling  oceanic  lithosphere  accretion.  For  example,  it  has 
been  shown  that  gravity  and  seafloor  depth  vary  systematically 
along  individual  spreading  segments  [e.g.,  Kuo  and  Forsyth, 
1988;  Prince  and  Forsyth,  1988;  Lin  et  ai,  1990;  Detrick  et  ai, 
1995]  and  appear  to  be  spreading-rate-dependent  [Parmentier 
and  Phipps  Morgan,  1990;  Un  and  Phipps  Morgan,  1992;  Sparks 
et  al.,  1993].  Such  variations  in  gravity  and  bathymetry  indicate 
segment-scale  changes  in  crustal  thickness  and/or  mantle  density 
and  thus  may  reflect  anomalies  in  along-axis  mantle 
temperatures.  Near  hot  spots,  however,  the  extent  of  along-axis 
variation  in  density  structure  is  broader  than  individual  ridge 
segments,  indicating  a  larger  scale  influence  by  mantle  plumes 
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[Anderson  et  al.,  1973;  Cochran  and  Talwani,  1977;  Bell  and 
Buck.  1992].  The  influence  of  mantle  plumes  on  crustal 
composition  is  also  evident  by  enrichments  of  trace  elements  and 
isotopes  along  the  Reykjanes  Ridge  near  the  Iceland  hot  spot, 
[Hartetal,  1973;  5c/iz7/mg,  1973,  1975a;  Vink,  1 984],  along  the 
Mid-Atlantic  Ridge  near  the  Azores  hot  spot  [White  et  al.,  1975, 
1976;  Schilling,  1975b],  and  along  the  Galdpagos  spreading 
center  near  the  Galapagos  hot  spot  [Schilling  et  ai,  1976,  1982; 
Verma  and  Schilling,  1982;  Verma  et  ai,  1983]. 

The  Galapagos  spreading  center  is  an  excellent  example  of  an 
oceanic  ridge  influenced  by  a  nearby  hot  spot.  At  present-day, 
the  spreading  center  lies  -170  km  north  of  the  Galapagos  hot  spot 
and  separates  the  Cocos  Plate  to  the  north  and  the  Nazea  Plate  to 
the  south  with  a  full  spreading  rate  of  4.5-6.8  cm/yr  [DeMets  et 
ai.  1990]  (Plate  la).  Spreading  segments  of  the  Galapagos 
spreading  center  trend  east-west  and  are  adjoined  by  north-south 
trending  transform  faults.  Hey  [1977]  proposed  that  the 
Galapagos  hot  spot  began  forming  the  Cocos  and  Carnegie 
Ridges  -20  Ma  and  then  migrated  southwest  with  respect  to  the 
Cocos  Plate  as  it  9ontinued  accreting  the  Cocos  Ridge.  The 
spreading  center  crossed  over  the  hot  spot  5-10  Ma  as  the 
Galapagos  Archipelago  began  its  formation  on  the  Nazea  Plate. 

As  for  the  present-day  interaction  between  the  hot  spot  plume 
and  spreading  center,  it  was  first  suggested  by  Morgan  [1978] 
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Plate  1.  (a)  Tectonic  map  of  the  Galapagos  spreading  system  encompassing  the  study  region  (rectangular  box). 
The  solid  dark  lines  mark  the  ridge  axis,  and  the  arrows  show  the  estimated  absolute  plate  motion  relative  to  the  hot 
spot  reference  frame,  (b)  Color-shaded  map  of  shipboard  and  DBDB5  bathymetiy  illuminated  from  the  north  and 
contoured  at  500-m  intervals.  Depths  shallower  than  1.6  km  arc  colored  red,  while  those  deeper  than  3.6  km  arc 
colored  violet.  Grid  spacing  is  5-min.  The  .spreading  center  is  marked  by  solid  white  lines  and  the  gravity  ship 
tracks  arc  marked  by  white  dotted  lines,  (c)  Color  map  of  free- air  gravity  along  ship  tracks  with  contour  interval  of 
10  mGal  and  gridded  with  5-min  spacing.  Gravity  values  >20  mGal  are  colored  red,  while  those  <-30  mGal  arc 
colored  violet.  The  contours  arc  drawn  from  interpolated  values  between  ship  tracks  and  are  masked  in  regions 
with  no  data. 
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that  Galiipagos  plume  material  feeds  through  a  mantle  conduit 
into  the  Galdpagos  spreading  center  giving  rise  to  the  Wolf- 
Darwin  seamount  chain  (Plate  lb).  Plume-fed  mantle  flow  along 
axis  was  suggested  by  Vogt  [1976]  to  explain  the  uniform 
increase  in  along-axis  topography  toward  the  hot  spot.  Further 
evidence  for  plume  flow  toward  and  along  the  spreading  center  is 
rare  earth  enrichments  along  the  ridge,  first  documented  by 
Schilling  et  ai  [1976].  Subsequent  studies  of  along-axis  variation 
in  rare  earth  element  and  isotopic  ratios  support  ideas  of  plume 
entrainment  to  the  ridge  axis  and  along-axis  dispersion  of  plume 
material  [Verma  and  Schillings  1982;  Verma  et  ai,  1983; 
Schillings  1985]. 

In  this  paper  we  present  evidence  for  a  regional  mantle 
temperature  anomaly  and  an  associated  crustal  thickness  variation 
beneath  the  Galapagos  spreading  center  imposed  by  the 
Galapagos  hot  spot.  We  first  isolate  the  subsurface  component  of 
the  gravity  field  by  making  topographic  and  crust-mantle 
interface  corrections  (the  mantle  Bouguer  corrections).  We  next 
examine  topographic  compensation  mechanisms  both  on-  and 
off-axis  by  comparing  theoretical  and  observed  mantle  Bouguer 
gravity  anomalies  along  isochrons  for  models  of  compensation 
from  crustal  thickness  variations  (Airy  compensation)  and 
compensation  from  laterally  varying  mantle  densities  (Pratt 
compensation).  Given  the  constraints  on  the  depth  of 
compensation,  we  then  examine  models  of  crustal  and  mantle 
density  structure  to  constrain  mantle  temperatures  along  the 
present-day  Galapagos  spreading  center.  Finally,  we  discuss  the 
temporal  evolution  of  axial  mantle  temperatures  in  the  past  7.7 
m.y.  and  its  implications  for  the  evolution  of  this  hot  spot-ridge 
system. 

Data 

Our  approach  for  investigating  mantle  temperature  anomalies 
requires  accurate  constraints  on  subsurface  density  structure 
which  is  reflected  directly  by  gravity  and  seafloor  topography. 
The  gravity  and  bathymetry  data  we  use  are  shipboard  data 
obtained  from  the  National  Geophysical  Data  Center  and  the 
Lamont-Doherty  Earth  Observatory.  The  data  set  covers  a  region 
within  84.4®-98.1®W  and  3.0°S-4.5°N,  encompassing  the 
Galapagos  spreading  center  and  the  Galapagos  hot  spot  (Plate 
lb).  We  also  use  high-resolution  gravity  and  bathymetry  data 
from  a  dense  survey  around  the  95.5°W  propagating  rift  tip 
[Phipps  Morgan  and  KleinrocK  1991].  The  bathymetry  data  are 
shipboard  depth  soundings  supplemented  with  digital  bathymetry 
(DBDB5)  between  ship  tracks.  DBDB5  data  points  within  5  min 
of  a  ship  data  point  are  eliminated  and  the  combined  data  set  is 
regridded  with  5-min  grid  spacing  to  produce  the  bathymetry  map 
shown  in  Plate  lb.  A  regional  bathymetric  swell  encompassing 
the  Galapagos  Archipelago  spans  -1300  km  along  the  ridge  axis. 
The  swell  shallows  along  the  ridge  axis  toward  91°W  by  1.1  km 
and  peaks  near  the  Galdpagos  hot  spot  which  is  centered  beneath 
the  island  Femandina  [White  et  ai.  1993}  (see  Figure  lb  for 
along-axis  profile). 

In  order  to  improve  the  internal  consistency  of  the  gravity  data 
we  use  the  method  of  Prince  and  Forsyth  [1984]  to  minimize 
discrepancies  in  gravity  measurements  at  ship  track  crossings. 
Applying  the  appropriate  DC  shifts  to  straight-line  track  segments 
reduces  the  total  rms  crossover  error  from  1 1 .2  to  5.5  mGal.  The 
value  of  5.5  mGal  is  therefore  our  estimate  for  data  uncertainty. 
After  applying  these  corrections  we  produce  the  5-min  grid  of 
free-air  gravity  shown  in  Plate  Ic.  In  this  map  we  observe  short- 
wavelength  (<100  km)  peaks  coinciding  with  topographic  highs; 
the  lowest  free- air  gravity  (-90  mGal)  occurs  over  the  flexural 


moat  of  the  Galapagos  Archipelago  and  the  highest  (+60  mGal) 
occurs  over  the  southeastern  flank  of  the  island  of  Isabela  (see 
Plate  lb  for  location  of  Isabela).  The  total  variation  in  free-air 
gravity  along  the  ridge  axis  is  -40  mGal. 

We  use  only  shipboard  gravity  rather  than  satellite- altimetry- 
derived  gravity  because  the  released  satellite  data  coverage  in  this 
region  is  still  sparse  and  the  shipboard  gravity  is  more  accurate. 
The  other  reason  for  using  only  shipboard  gravity  concerns  the 
accuracy  of  topographic  corrections  which  rely  on  accurate 
bathymetric  measurements.  Since  shipboard  gravity  and 
bathymetric  measurements  are  taken  at  the  same  points, 
topographic  corrections  to  the  free-air  gravity  are  the  most 
accurate  possible. 

Gravity  Data  Reduction 

A  significant  portion  of  the  free-air  gravity  is  caused  by 
seafloor  topography.  Therefore,  to  investigate  subsurface  density 
structure,  to  which  we  will  relate  mantle  temperature  anomalies, 
we  apply  a  mantle  Bouguer  correction.  Using  Parker’s  [1973] 
spectral  method,  we  subtract  from  the  free-air  gravity  the 
theoretical  gravity  signal  of  the  seafloor-water  interface  and 
crust-mantle  (Moho)  interface  assuming  a  crustal  layer  of 
constant  thickness  (6  km)  and  density  (2800  kg/m^).  We  take  the 
density  of  the  mantle  to  be  3300  kg/m^. 

The  resulting  mantle  Bouguer  anomaly  (MBA)  shows  that 
most  of  the  short- wavelength  (<100  km)  variations  caused  by 
local  topography  are  removed,  leaving  a  broad  oval-shaped 
negative  anomaly  aligned  along  the  spreading  axis  between 
~97°W  and  ~85°W  (outlined  by  yellow  shading,  Plate  2a). 
Superimposed  on  this  oval-shaped  anomaly  are  high-amplitude 
negative  branches  over  the  Cocos  Ridge  (<-100  mGal)  and 
Galdpagos  Archipelago  (<-300  mGal)  reflecting  the  thickened 
crust  of  these  edifices.  Along  the  ridge  axis,  10-20  mGal 
variations  in  MBA  occur  at  segmentation  length  scales  (100-200 
km),  but  the  most  prominent  feature  is  the  long-wavelength 
decrease  by  ~90  mGal  along  axis  toward  91®W  (Plate  2b).  For 
comparison  with  other  oceanic  spreading  centers,  this  90-mGal 
anomaly  is  approximately  twice  the  segmentation-scale  MBA 
variation  along  the  Mid- Atlantic  Ridge  [e.g.,  Lin  et  ai,  1990; 
Detrick  et  ai.  1995]  and  about  10  times  the  MBA  variation  along 
the  East  Pacific  Rise  (8.8°-13.5°N)  [Madsen  et  ai,  1990]. 

The  minimum  in  MBA  occurs  near  9l°W  on  the  southern 
segment  of  the  91°W-transform  offset,  which  is  also  the  point  of 
the  ridge  axis  closest  to  the  Galapagos  hot  spot  (point  P,  Plate 
2a).  The  decrease  in  MBA  is  nearly  symmetric  about  point  P  and 
uniform  along  the  650-km  ridge  length  to  each  side  of  point  P. 
This  wavelength  is  comparable  in  extent  to  topographic  swells  of 
other  hot  spots  such  as  Hawaii  (-1500  km  across  the  island 
chain).  Cape  Verde  (-1500  km  in  diameter)  [Crough.  1983],  and 
Iceland  (-2000  km  in  diameter)  [White,  1988]. 

Comparison  of  this  along-axis  MBA  with  along-axis  variations 
in  bathymetry  and  basalt  chemistry  reveals  a  close  correlation 
between  the  four  anomalies  (Figure  1).  All  anomalies  peak  at  or 
near  point  P,  all  extend  over  comparable  length  scales,  and  all 
decrease  in  amplitude  nearly  symmetrically  eastward  and 
westward  away  from  point  P.  The  peak  in  the  La/Sm  anomaly 
coincides  with  that  of  K2O.  MgO,  and  a  minimum  in  FeO 
[Schilling  et  ai,  1982],  while  the  peak  in  *‘'Sr/®^Sr  coincides  with 
a  minimum  in  '“^^Nd/'^^Nd  [Verma  et  ai,  1983].  Both 
geochemical  signatures  are  attributed  to  a  source  heterogeneity 
associated  with  the  Galdpagos  hot  spot  [Verma  et  ai,  1983]. 
Although  the  peak  in  “^Sr/^'^Sr  occurs  -100  km  west  of  point  P, 
this  offset  is  small  relative  to  the  total  wavelength  of  the  above 
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Figure  1.  (a)  Along-axis  mantle  Bouguer  (MBA)  and  bathymetry  profiles  are  compared  with  along-axis  variations 
(b)  in  [La/Sm]^fand  *"^Sr/®^Sr  \Verma  et  al,  1983].  Note  that  the  peaks  for  all  anomalies  except  for  ^"^Sr/^^Sr 
coincide  at  point  P. 


anomalies  and  within  the  150  to  300-km  diameter  suggested  for 
intraplate  hot  spots  [Epp,  1984;  McNutt,  1989].  The  correlation 
of  MBA  and  bathymetry  with  basalt  chemistry  suggests  that  the 
along-axis  density  structure  is  closely  related  to  the  enriched 
material  introduced  by  Galapagos  plume  beneath  the  ridge  axis. 

The  final  step  in  our  gravity  analysis  is  to  remove  the 
predictable  effects  due  to  lithospheric  cooling.  Calculation  of  the 
three-dimensional  (3-D)  distribution  of  temperature-dependent 
mantle  densities  for  passive  mantle  upwelling  is  relatively  simple 
using  a  standard  method  first  presented  by  Phipps  Morgan  and 
Forsyth  [1988].  We  use  a  spectral  method  to  solve  for  flow  of  a 
constant-viscosity  mantle,  driven  by  two  spreading  plates  with 
the  geometry  of  the  Galapagos  spreading  center.  Using  finite 
difference  approximations  for  the  conductive-advective  heat 


balance  equation,  we  solve  for  steady  state  mantle  temperatures, 
from  which  we  derive  mantle  densities  assuming  a  thermal 
expansion  coefficient  of  3.4x10'^  The  integrated  gravity 
fields  from  each  density  layer  down  to  a  100-km  depth  represent 
the  lithospheric  cooling  contribution  to  the  gravity  field  [Kuo  and 
Forsyth,  1988]. 

We  subtract  the  lithospheric  cooling  effects  from  the  MBA  to 
produce  the  residual  mantle  Bouguer  gravity  anomaly  (RMBA), 
shown  in  Plate  2c.  The  oval-shaped  MBA  low,  observed  between 
~9TW  and  -85®W,  is  removed  by  the  lithospheric  cooling 
correction;  what  remain  are  high-amplitude  negative  anomalies 
branching  from  the  ridge  axis,  over  the  Galapagos  Archipelago 
and  the  Cocos  and  Carnegie  Ridges.  These  negative  RMBA 
branches  reflect  the  anomalous  volcanic  and  mantle  density 


Plate  2.  (a)  Contour  and  color-shaded  map  of  the  mantle  Bouguer  anomaly.  Anomaly  values  >20  mGal  are  colored 
red,  while  those  <-120  mGal  are  colored  violet.  Interpolated  values  between  ship  tracks  are  masked,  and  the 
spreading  center  and  islands  are  marked  in  white.  Note  the  oval-shaped  negative  anomaly  aligned  along  the 
spreading  center  between  ~97°W  and  ~85°W  (outlined  by  the  yellow  shading)  and  the  negative  anomaly  branches 
of  the  Cocos  Ridge  and  Galapagos  Archipelago.  The  five  white  profiles  north  of  the  spreading  center  mark 
isochrons  from  Wilson  and  Hey  [1995]  used  for  the  off-axis  analysis.  Profiles  are  dashed  in  regions  where 
magnetic  lineations  are  extrapolated,  (b)  Mantle  Bouguer  gravity  values  extracted  along  the  spreading  center. 
Note  that  the  anomalies  reach  a  minimum  at  point  P.  where  the  ridge  axis  is  closest  to  the  hot  spot.  The  arrows 
mark  locations  of  transform  offsets,  (c)  Map  of  residual  mantle  Bouguer  anomaly  with  contour  interval  of  20  mGal 
and  a  color  range  of  -90  to  +50  mGal.  Note  high-amplitude  negative  anomalies  along  the  Cocos  Ridge,  the 
Darwin-Wolf  seamounts,  and  the  Galapagos  Islands  (shown  in  white),  (d)  Along-axis  profile  of  residual  mantle 
Bouguer  anomaly  showing  -100  mGal  decrease  from  the  east  and  west  toward  point  P.  Arrows  mark  transform 
offsets. 
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structure  left  along  the  Galapagos  hot  spot  tracks.  The  dominant 
effect  of  the  lithospheric  correction  on  the  along-axis  profiles  is 
to  reduce  the  segment-scale  variations  in  the  MBA;  it  does  not 
appreciably  affect  the  long- wavelength  decrease  due  to  the  hot 
spot  (Plate  2d).  Although  the  amplitude  of  the  along-axis  RMBA 
is  increased  slightly  (by  10  mGal)  from  that  of  the  along-axis 
MBA,  the  lateral  extent  and  location  of  the  maximum  are  the 
same  for  the  two  anomalies.  For  this  reason,  we  focus  on  the 
MBA  for  the  remainder  of  the  paper. 

Compensation  of  Topography 

The  mantle  Bouguer  correction  has  been  widely  used  as  a  first- 
order  correction  for  oceanic  crustal  structure  [e.g.  Kuo  and 
Forsyth^  1988;  Lin  et  al,  1990;  Madsen  et  al,  1990;  Blackman 
and  Forsyth,  1991;  Morris  and  Detrick,  1991]  since  the  total 
global  variation  in  oceanic  crustal  thickness  is  ~3  km  about  a 
mean  of  6  km  [Chen,  1992].  Our  assumption  of  a  constant  6-km- 
thick  crust  is  merely  a  first  approximation  of  crustal  structure 
from  which  we  reference  departures  in  density  structure. 
Deviations  from  this  reference  model  could  be  lateral  variations 
in  crustal  thickness,  lateral  mantle  density  variations,  or  a 
combination  of  the  two.  The  nonuniqueness  of  gravity  solutions 
necessitates  additional  constraints.  An  obvious  constraint  is 
topography  since,  if  in  isostatic  equilibrium,  it  too  depends 
directly  on  density  structure.  Here  we  test  two  modes  of  isostatic 
compensation:  (1)  crustal  compensation  (Airy  isostasy)  and  (2) 
compensation  from  lateral  density  variations  in  the  mantle  (Pratt 
isostasy). 

Airy  and  Pratt  Isostasy  Admittance  Functions 

The  theoretical  mantle  Bouguer  gravity  anomaly  due  to  the 
two  modes  of  isostatic  compensation  is  solved  using  standard 
spectral  methods  as  follows.  In  the  spectral  domain,  the  mantle 
Bouguer  gravity  anomaly  5(k)  is  related  to  bathymetry  HQi)  by 
an  isostatic  response  function,  or  admittance  function  Q{k), 
according  to 

B(k)=e(A)-/f(k)  (I) 

where  k  is  the  two-dimensional  wavenumber,  /t=lkI=2rc/X. 
Included  in  Q{k),  are  the  effects  of  density  structure  that  differs 
from  the  “reference”  structure  (i.e.,  a  crust  of  uniform  thickness 
overlying  a  mantle  of  uniform  density).  In  Airy  compensation 
models,  it  is  the  crustal  structure  that  differs  from  the  “reference” 
since  topography  is  assumed  to  be  supported  by  laterally  varying 
crustal  thickness.  If  we  assume  that  elevated  topography  is 
supported  by  crust  that  is  anomalously  thick,  the  admittance 
function  must  include  two  terms  to  account  for  the  effects  at  two 
interfaces  as  follows: 

l2(^)=-27cG[Apexp(-^2^)+p^exp(-A2^)],  (2) 

where  (7=6.67x10’^^  m^/kg  is  Newton^s  gravitational  constant, 
Ap  is  the  crust-mantle  density  contrast  (500  kg/m^),  and  p^  is  the 
crust- water  density  contrast  (1800  kg/m^).  The  first  term  replaces 
the  attraction  of  mantle  by  that  of  crast  at  the  “reference”  Moho 
(Zc=8.7  is  the  average  seafloor  depth  of  2.7  km  plus  6.0  km).  The 
second  term  accounts  for  the  effects  of  the  Airy  crustal  root  at  its 
assumed  mean  depth  of  1 1.0  km  beneath  the  sea  surface. 

For  Pratt  compensation,  it  is  mantle  density  that  differs  from 
the  “reference”  structure  since  topography  is  assumed  to  be 
compensated  by  laterally  varying  mantle  densities.  The 


amplitude  of  the  density  reduction  in  a  vertical  column  required 
to  support  a  given  topographic  elevation  is  HpJZp  where  p^  is 
the  mantle-water  density  contrast  (2300  kg/m^)  and  Zp  is  the  depth 
of  compensation.  The  Pratt  admittance  function  must  therefore 
consider  the  integrated  effects  of  all  density  layers  from  to  (z^+ 
Zp)  and  is  thus  defined  as 

Q{k)-  -27tCp^  ^[1  exp(  kZp)] 

kzp 

Results 

Mantle  Bouguer  profiles  taken  along  the  present-day  ridge  axis 
and  selected  isochrons  (Plate  2a)  are  compared  with  the  Airy  and 
Pratt  theoretical  profiles  (Figure  2a).  The  standard  deviation 
misfit  between  theoretical  and  observed  profiles  is  plotted  versus 


Age  (Ma) 

Figure  2.  (a)  Profiles  of  observed  mantle  Bouguer  anomalies 
(shaded  profiles)  are  compared  with  theoretical  models  for 
different  assumed  compensation  mechanisms  and  depths.  The 
locations  of  the  off-axis  crustal  isochrons  (labeled  with  ages  from 
Wilson  [1993])  are  in  Plate  2c.  (b)  The  standard  deviation  misfit 
is  plotted  against  crustal  age  for  the  five  compensation  models 
tested.  Note  that  models  of  shallow  depths  of  compensation  yield 
the  smaller  misfits. 
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age  in  Figure  2b.  For  the  crustal  ages  examined  (0.0-7.67  Ma) 
the  standard  deviation  misfit  for  the  Airy  compensation  profiles 
has  a  nearly  constant  value  of  -6  mGal,  which  is  close  to  our 
estimated  data  uncertainty  of  5.5  mGal.  For  the  Pratt 
compensation  models,  standard  deviation  misfits  increase  with 
compensation  depth  and  with  age.  Along  the  present-day 
axis,  the  Pratt  models  of  z^=50  and  100  km  are  the  most 
reasonable  with  standard  deviations  of  5.9  and  7.1  mGal, 
respectively. 

Although  the  Airy  profiles  yield  the  lowest  misfits  to  the 
observed  MBA,  the  Pratt  calculations  with  shallow  compensation 
depths  (Zp=50  and  100  km)  also  yield  small  misfits.  Most  of  the 
misfit  by  the  shallow  Pratt  calculations  appears  to  be  due  to  short- 
wavelength  variations  (<200  km)  in  the  observed  MBA  which 
may  come  from  local  variations  in  crustal  structure.  We  thus  do 
not  exclude  the  possibility  that  at  least  some  of  the  gravity  and 
topographic  signal  is  due  to  density  variations  in  the  shallow 
mantle.  The  increase  in  misfits  with  age  for  the  Pratt  calculations 
may,  however,  reflect  a  decrease  in  the  mantle  contribution 
relative  to  that  of  the  crust  along  paleoridge  axes. 

Topography  may  also  be  supported  dynamically  by 
lithospheric  or  upper  mantle  stresses.  Previous  work  has  shown 
that  shallow  stresses  induced  by  plate  spreading  contribute 
significantly  to  axial  topography  [e.g.,  Phipps  Morgan  et  ai, 
1987;  Lin  and  Parmentier,  1989;  Small  and  Sandwell,  1989; 
Chen  and  Morgan,  1990;  Neumann  and  Forsyth,  1993]. 
Neumann  and  Forsyth  11993],  for  example,  demonstrated  that  the 
correlation  between  gravity  and  topography  along  the  Mid- 
Atlantic  Ridge  is  due  to  dynamic  stresses  in  the  lithosphere  which 
depend  on  crustal  thickness  and  mantle  thermal  structure. 
However,  for  the  Gal^ipagos  spreading  rates  of  48-64  mm/yr  and 
possible  crustal  thickness  structure,  these  extension-related 
stresses  would  be  small  [Neumann  and  Forsyth,  1993]. 
Significant  topography  (>1  km)  can  also  be  maintained  by 
viscous  stresses  in  a  convecting  mantle  [Anderson  et  ai,  1973; 
McKenzie  et  ai,  1980].  If  viscous  stresses  are  important  along 
the  Galapagos  spreading  center,  they  must  be  associated  with  low 
densities  in  the  shallow  mantle  as  indicated  by  our  MBA 
modeling  (Figtire  2).  We  thus  suggest  that  the  long-wavelength 
topography  of  the  Galapagos  spreading  center  and  nearby  Cocos 
Plate  is  essentially  isostatic  and  is  supported  by  density  anomalies 
primarily  within  100  km  beneath  the  seafloor. 

Present-Day  Axial  Mantle  Temperatures 

As  demonstrated  above,  crustal  thickness  and  shallow  mantle 
density  variations  are  both  likely  sources  of  topographic 
compensation.  We  suggest  that  they  both  occur  and  that  both  are 
controlled  by  mantle  temperature:  crustal  thickness  by 
temperature-enhanced  melting,  and  mantle  density  by  thermal 
expansion.  For  the  following  analysis,  we  investigate  the  mantle 
temperature  variation  required  to  generate  the  ~90-mGai 
variation  in  along-axis  MBA  and  the  -1.1 -km  increase  in  axial 
topography. 

Model  Configuration 

Using  the  same  numerical  methods  as  were  used  for  the 
lithospheric  cooling  calculations,  we  solve  for  3-D  mantle 
temperatures  due  to  passive  upwelling;  this  time,  however,  we 
impose  a  temperature  anomaly  AT  at  the  base  layer  (Figure  3). 
To  estimate  the  additional  crust  that  may  result  from  a  given  AT, 
we  take  the  fraction  of  partial  melting  /  to  depend  on  mantle 
temperature  rbyy=(T-T,)/6()0®,  where  T,  is  the  mantle  solidus 


Longitude 


Figure  3.  Schematic  illustration  of  our  simplified  3-D  melt 
generation  models.  The  ridge  axis  (shaded  gray  lines  at  the  top 
surface)  is  offset  by  the  91®  transform  fault.  The  region  of  melt 
transport,  with  width  25,  is  outlined  by  the  bold  dark  lines.  The 
melting  zone  beneath  ridge  segments  is  shown  on  depth  cross 
sections  as  shaded  triangular  shapes;  curved  arrows  denote  melt 
transport  toward  the  ridge  axis.  The  region  of  melt  transport  at 
the  base  layer  is  shaded  gray  according  to  the  imposed 
temperature  with  the  greatest  temperature  anomaly  near  the  91® 
transform  fault. 

and  600°C  is  the  supersolidus  temperature  required  to  fully  melt  a 
unit  volume  of  peridotite  [Reid  and  Jackson,  1981].  The  rate  of 
melt  generation  /  depends  on  the  gradient  of/ and  mantle  flow 
rate  v  by  /  =v*V/  [Reid  and  Jackson,  1981].  We  estimate  the 
mantle  solidus  to  be  linearly  dependent  on  pressure  and  thus 
depth  z  (in  kilometers  below  the  seafloor),  by  7^=1 100®C  + 
3.25(®C/km)z.  Crustal  accretion  at  the  ridge  crest  depends  on  the 
spatial  distribution  of  melting  and  subsequent  migration  of  melt 
toward  the  ridge.  The  process  of  melt  migration  is  still  largely 
unconstrained;  therefore  we  simplify  this  calculation  by  treating 
ridge  segments  as  line  sinks  which  draw  in  melt  from  the  mantle 
below  [Phipps  Morgan  and  Forsyth,  1988].  Assuming  that  melt 
migrates  over  a  finite  lateral  extent,  we  define  a  width  5,  on  each 
side  of  the  ridge  axis,  as  the  region  of  melt  transport  (Figure  3); 
outside  of  this  melt  transportation  zone,  we  assume  that  melts  are 
carried  away  by  the  cooling  lithospheric  plates  thus  do  not 
contribute  to  the  crust.  We  also  assume  that  a  small  melt  fraction 
ffy  is  retained  in  the  mantle  matrix  within  this  zone  of  melt 
transport.  We  adjust  the  parameters  5  and  /<,  such  that  the 
resulting  crustal  thickness  for  a  normal  base  layer  temperature 
(1350®C)  is  6  km  at  the  ridge  segment  centers.  This  result  is 
achieved  for/,  values  of  0-6%  [Forsyth,  1992]  and  corresponding 
5  values  of  30-50  km.  We  assume  that/^  does  not  vary  along- 
axis  therefore  it  does  not  contribute  to  the  long-wavelength 
variation  in  mantle  density.  The  base  layer  is  set  to  a  depth  of 
160  km  to  ensure  that  the  entire  melting  region  is  included. 

Base  Layer  Temperature  Distributions 

For  the  base  layer  temperature  anomalies,  we  investigate  three 
geometries.  In  our  first  set  of  calculations  (model  A),  we  vary 
temperatures  linearly  along-axis  with  the  maximum  anomaly 
beneath  the  91°W  transform  (Figure  4a).  This  is  the  simplest 
model,  designed  to  test  the  effects  of  strictly  along-axis  variation 
in  temperatures.  For  the  second  set  of  calculations  (model  B),  we 
impose  a  circular  anomaly  centered  on  the  island  Femandina, 
thought  to  mark  the  current  location  of  the  Galapagos  hot  spot 
center  [White  et  ai,  1993]  (point  H,  Figure  4b).  Temperatures 
decrease  linearly  away  from  point  H.  We  envision  this  model  to 
represent  the  temperature  distribution  from  a  radial  dispersion  of 
plume  material  from  the  hot  spot  center.  For  the  third  set  of 
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Figure  4.  (a)  Map  of  temperatures  imposed  at  the  base  layer  of  our  best  fitting  linear  model  A.  Temperatures  in 
the  zone  of  melt  transport  are  shaded  to  emphasize  the  importance  of  this  region  in  gravity  and  bathymetry 
calculations.  Point  H  marks  the  location  of  the  Galapagos  hot  spot,  while  point  P  marks  the  location  of  the  peak  in 
along-axis  MBA  and  baihymetiy.  (b)  Temperatures  imposed  at  the  base  layer  of  our  best  fitting  circular  model 
(model  B)  with  maximum  temperature  located  at  point  H.  Arrows  denote  hypothetical  radial  dispersion  of  hot  spot 
material  from  point  H.  (c)  Base  layer  temperatures  imposed  for  our  best  fitting  elliptical  model  (model  C).  Arrows 
denote  plume  channeling  from  point  H  to  point  P  (arrow  1 )  and  then  along-axis  away  from  P  (arrows  2). 


calculations  (model  C),  we  use  an  elliptical  temperature  anomaly 
which  is  centered  midway  along  the  9rw  transform  and 
decreases  linearly  away  from  the  ellipse  center  (Figure  4c).  This 
model  is  designed  to  approximate  the  temperature  distribution 
along  the  ridge  axis  that  might  result  from  Schilling's  [1985, 
1991]  plume  flow  model  which  incorporates  ideas  of  Vogt  [1976] 
and  Morgan  [1978].  According  to  this  model,  Galapagos  plume 
material  feeds  through  a  conduit  connecting  point  H  to  the  ridge 
axis  (arrow  1 ,  Figure  4c),  and  then  disperses  preferentially  along 
axis  (2  arrows.  Figure  4c).  We  approximate  the  preferential 
along-axis  flow  as  an  ellipse  aligned  with  the  ridge  axis.  Each  of 
the  three  models  has  a  base  layer  temperature  maximum  near 
point  P  directly  beneath  the  ridge  axis  with  a  gradual  decrease 
along  axis  toward  the  east  and  west  edges  of  the  study  region. 
Gravity  and  bathymetry  calculations  for  these  models  are 
sensitive  mostly  to  temperature  conditions  within  the  region  of 
melt  transport  since  only  melts  in  this  region  are  assumed  to 
contribute  to  accretion  of  the  crust. 

Results 

Gravity  calculations  for  the  three  models  are  done  by  applying 
Parker's  [1973]  method  to  the  density  layers  in  the  mantle  and 
the  crust-mantle  interface  treating  the  crustal  thickness  as  only 
varying  along-axis.  Since  we  have  shown  that  the  long- 
wavelength  seafloor  topography  is  compensated  at  shallow 


depths,  we  calculate  theoretical  bathymetry  assuming  Airy 
compensation  for  the  crust  and  Pratt  compensation  for  the  mantle 
shallower  than  160  km.  Figures  5a  and  5b  compare  theoretical 
results  of  model  A  for  different  base  layer  temperature  anomalies  ‘ 
at  point  P  (A7^)  with  observed  MBA  and  bathymetry  profiles. 
Profile  sections  near  transform  faults  are  omitted  since  our 
models  gave  unrealistipally  small  crustal  thicknesses  due  to  local 
cooling  effects  near  ridge  segment  ends.  Removal  of  these  local 
effects,  however,  do  not  affect  the  larger-scale  thermal  anomaly 
related  to  the  Galapagos  hot  spot. 

As  illustrated  in  Figures  5a  and  5b,  the  model  for  AT],  of  50°C 
best  fits  both  the  MBA  and  bathymetry  profiles.  The  A7],=25  and 
75°C  solutions  are  the  upper  and  lower  limits  for  model  A.  Table 
1  outlines  the  corresponding  results  of  models  B  and  C  and  the 
associated  standard  deviation  misfits.  Despite  differences  in  the 
detail  3-D  temperature  structure  between  the  three  models,  all 
three  suggest  similar  values  of  AT^  with  comparable  minimum 
misfit.  This  finding  indicates  that  axial  crustal  and  density 
structure  is  sensitive  primarily  to  temperatures  directly  beneath 
the  ridge  axis  and  insensitive  to  subtle  temperature  changes  away 
from  the  ridge  axis.  We  conclude  AT^  to  be  ~50±25°C  with  a 
corresponding  crustal  thickening  of  3±1  km  (Figure  5c).  As  the 
crust  thickens  toward  point  P,  it  gives  rise  to  70-75%  of  the 
topographic  swell  and  MBA  gravity  signal.  The  remaining  25- 
30%  of  topography  and  gravity  signal  is  supplied  by  the 
anomalously  hot  and  less  dense  mantle  beneath  the  ridge  axis. 
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Figure  5.  The  observed  profiles  (small  crosses)  of  (a)  along-axis  MBA  and  (b)  bathymetry  are  compared  with 
theoretical  profiles  from  model  A  with  different  values  of  point  P,  base  layer  temperature  anomalies,  AT^.  Sections 
of  profiles  near  transform  faults  are  omitted  to  accentuate  the  broad  wavelength  anomaly  (solid  lines)  associated 
with  the  hot  spot.  The  best  fitting  profiles  are  denoted  by  bold  lines,  (c)  Predicted  along-axis  cross  section  of  the 
crustal  structure  based  on  model  A  results.  The  Moho  is  drawn  according  to  our  best  fitting  result  (+50'^C  model); 
the  two  other  profiles  are  drawn  according  to  the  +25°C  (shallower  curve)  +75®C  (deeper  curve)  results.  The 
Moho  boundary  is  omitted  near  transform  faults  as  marked  by  arrows.  Densities  are  in  grams  per  cubic  centimeter. 


Our  crustal  model  is  consistent  with  estimates  of  Feighner  and 
Richards  [1994]  based  on  flexural  modeling  of  gravity  near  the 
Galapagos  Archipelago.  Confirming  this  crustal  model,  however, 
requires  future  marine  seismic  experiments  since  few  seismic 
constraints  exist  to  date. 

Discussion 

Our  primary  focus  in  this  study  is  the  effects  of  mantle 
temperature  on  crustal  thickness  and  on  mantle  density  changes 
by  thermal  expansion.  A  number  of  factors  not  incorporated  into 
our  models  may  also  contribute  to  crustal  thickness  and  mantle 
density  structure  and  lead  to  changes  in  our  AT^  estimate.  These 
include  (1)  compositionally  dependent  and  disequilibrium 
melting,  (2)  melt  depletion  and  latent  heat  loss  in  the  mantle,  (3) 
buoyancy-driven  mantle  flow,  and  (4)  mantle  compositional 


effects  on  melting  and  on  mantle  densities.  We  briefly  discuss 
these  factors  below. 

Compositionally  dependent  and  disequilibrium  melting. 
The  simple  linear  melt  function  and  linear  depth-solidus  relation 
that  we  used  was  defined  by  Reid  and  Jackson  [1981]  based  on 
results  of  batch  melting  experiments  in  which  melt  maintained 
equilibrium  with  the  remaining  solid  phases.  If  melt  is  extracted 
rapidly  in  the  mantle  such  that  it  fails  to  equilibrate  with  the 
matrix,  then  the  solidus  of  the  depleted  residue  increases  with 
increasing  melt  extents  [Kinzler  and  Grove^  1992;  Cordery  and 
Phipps  Morgan,  1993].  If  this  disequilibrium  melting  scenario  is 
the  dominant  process  in  the  mantle,  then  a  greater  AT^  than  we 
estimated  may  be  required  beneath  the  Galapagos  spreading 
center  to  thicken  the  crust  by  3±1  km. 

Melt  depletion  and  latent  heat  loss  in  the  mantle.  In 
addition  to  inhibiting  melting,  melt  depletion  may  also  reduce 
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Table  1.  Model  Results 


MBA 

Bathvmetrv 

Best  Fit 
Ar^.-c 

Misfit, 

mGal 

Best  Fit 
AT,°C 

p 

Misfit, 

km 

Model  A 

50±25 

9 

50±25 

0.13 

Model  B 

49±25 

10 

49±25 

0.13 

Model  C 

47±25 

8 

49±25 

0.14 

mantle  densities,  primarily  by  decreasing  the  Fe/Mg  ratio  of  the 
residue  [Oxburgh  and  Parmentier,  1977].  The  opposite  effect, 
however,  may  result  from  latent  heat  removal  which  cools  the 
mantle  thus  increases  its  density.  Numerical  experiments  by 
Magde  et  al.  [1995]  indicate  that  to  generate  an  oceanic  crust  of 
normal  thickness,  the  two  factors  would  lead  to  a  net  decrease  in 
mantle  densities  of  the  order  of  1%,  or  ~6  times  the  thermal 
expansion  effect  of  heating  the  mantle  by  50°C.  This  potential 
density  reduction  may  contribute  to  the  Galapagos  bathymetry 
and  gravity  anomalies  significantly  enough  that  the  crustal 
thickness  and  thus  AT],  are  smaller  than  we  estimated. 

Buoyancy-driven  mantle  flow.  In  addition  to  their  direct 
signature  on  surface  observables,  mantle  density  variations  lead 
to  buoyancy  forces  which  drive  mantle  flow.  Beneath  normal 
oceanic  spreading  centers  the  dominant  sources  of  buoyancy  are 
most  likely  melt  depletion-related  and  melt  retention-related 
density  reductions  [Jha  et  al,  1994].  If  buoyancy  forces  are 
important,  they  are  likely  to  enhance  vertical  flow  and  increase 
the  rate  of  melting  and  thus  may  lead  to  a  lower  A7^  prediction. 

Mantle  compositional  effects  on  melting  and  on  mantle 
densities.  A  hot  spot-related  temperature  anomaly  as 
investigated  in  this  study  is  an  obvious  source  for  thickened  crust 
and  reduced  mantle  densities;  however,  mantle  source 
heterogeneity  may  also  play  important  roles.  Enrichment  in 
volatile  [Bonatti,  1990]  or  incompatible  elements  [Michael  et  al, 
1994]  in  the  mantle  may  enhance  melting  and  thus  yield  a  thicker 
crust  for  a  given  mantle  temperature  anomaly.  While  there  is 
little  evidence  for  a  volatile  enrichment  beneath  the  Galapagos 
spreading  center,  there  is  evidence  for  an  increase  in  incompatible 
element  concentration  toward  point  P  from  ridge  axis  samples 
[Schilling  et  al,  1982;  Langmuir  et  al,  1992].  In  addition,  a 
decrease  in  Fe/Mg  was  observed  in  axial  samples  toward  point  P 
[Schilling  etal,  1982;  Langmuir  et  al,  1992],  possibly  reflecting 
a  low  Fe/Mg  and  thus  low-density  mantle  source  near  the 
Galapagos  hot  spot.  Including  such  heterogeneities  of  the  mantle 
source  in  incompatible  element  content  and  Fe/Mg  ratio  may 
yield  a  lower  ATp  estimate. 

In  summary,  while  considering  factor  1  would  increase  an 
estimate  of  AT^,  considering  factors  2,  3,  and  4  would 
substantially  decrease  an  estimate  of  ATp.  We  therefore 
anticipate  that  our  estimate  of  ATp  is  an  upper  bound,  although 
the  interplay  of  the  above  four  factors  may  be  complex  and 
requires  comprehensive  investigation  that  is  beyond  the  scope  of 
this  study. 

By  considering  the  first-order  aspects  of  mantle  flow,  heat 
transport,  and  decompression  melting,  we  have  established  a 
relation  between  crustal  thickness,  temperature-dependent  mantle 
density,  and  mantle  temperature  anomaly.  Our  approach  is 
consistent  with  previous  studies  of  intraplate  hot  spot  anomalies. 


For  example,  Crough  [1978],  Sleep  [1990],  and  McNutt  [1987] 
constrained  hot  spot  temperature  anomalies  based  on  mantle- 
density  anomalies  which  they  took  to  be  strictly  temperature 
dependent.  McKenzie  [1984]  constrained  hot  spot  anomalies 
based  on  estimates  of  crustal  thickness  assuming,  as  we  do,  that 
melting  occurs  under  equilibrium  conditions.  Our  relationship 
between  AT  and  the  mantle  component  of  topography  is 
essentially  the  same  as  Sleep's  [1990],  and  our  relationship 
between  AT  and  crustal  thickness  is  consistent  with  that  of 
McKenzie  [1984]  (50-75°C  for  crustal  thickening  of  2-4  km). 

Our  constraints  on  AT  beneath  the  Galapagos  spreading  center 
have  implications  for  the  nature  of  heat  transport  both  along  the 
ridge  axis  and  from  the  hot  spot  to  the  ridge  axis.  Using  Feighner 
and  Richard’s  [1994]  estimate  for  the  volcanic  thickness  of  the 
Galapagos  Archipelago  (15-20  km)  and  McKenzie’s  [1984] 
melting  relationships,  we  estimate  a  temperature  anomaly  of 
~200°C  at  the  hot  spot  center  (point  H,  Figure  6).  This 
temperature  estimate  is  similar  to  the  214°C  anomaly  estimated 
by  Schilling  [1991]  based  on  buoyancy  flux  arguments. 
Considering  our  upper  (75°C)  and  lower  (25°C)  estimates  for  the 
temperature  anomaly  at  point  P,  the  average  gradient  from  the  hot 
spot  to  the  ridge  axis  (H  to  P,  Figure  6)  is  0.74- 1 .03°C/km.  In 
contrast,  the  along-axis  gradient  is  only  0.04-0.1  l®C/km. 
Therefore  any  successful  models  of  sublithospheric  plume 
dispersion  must  yield  an  along-axis  temperature  gradient  that  is 
an  order  of  magnitude  less  than  that  from  the  hot  spot  to  the  ridge 
axis.  Rigorous  investigation  of  this  question  requires  further 
experimental  [Kincaid,  1994]  and  numerical  [Rowley  et  al,  1992] 
work. 

Paleoaxial  Temperature  Anomalies 

To  better  constrain  the  temporal  evolution  of  the  Galapagos 
ridge-hot  spot  system,  we  next  examine  MBA  and  bathymetry 
anomalies  along  paleoaxes  of  the  Galapagos  spreading  center. 
From  our  model  calculations  we  derive  empirical  relations 
between  AT  and  MBA  and  bathymetry  that  we  then  use  to 
estimate  past  temperature  anomalies  from  the  observed 
amplitudes  of  along-isochron  MBA  and  bathymetry  anomalies. 
In  order  to  apply  relationships  derived  from  the  active  spreading 
center  to  off-axis  anomalies,  we  must  make  two  assumptions. 
First,  we  assume  that  any  off-axis  crustal  accretion  on  the  Cocos 
Plate  is  insignificant  and  that  the  spreading  rate  has  not  changed 
significantly  over  the  past  7.7  m.y.  Second,  we  assume  the 


Longitude 


Figure  6.  Map  of  the  Galapagos  region  marked  with  estimated 
base  layer  temperature  anomalies  at  various  points  along  the  ridge 
axis  (solid  line)  and  at  the  hot  spot  center  (point  H).  Arrows  point 
in  the  direction  of  decreasing  temperature  anomalies  from  the  hot 
spot  to  ridge-axis  (arrow  1)  and  along  the  ridge  axis  (arrows  2). 
Estimated  temperature  gradients  in  both  directions  are  labeled. 
Note  that  gradient  1  is  10-20  times  greater  than  gradient  2. 
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Figure  7.  The  empirical  relationships  between  base  layer  temperature  anomaly  and  (a)  MBA  and  (b)  bathymetry 
anomaly  (solid  lines).  Also  plotted  are  derived  AT  values  for  maximums  in  observed  MBA  and  bathymetry  along 
crustal  isochrons  (dots).  Errors  in  gravity  (12  mGal)  and  bathymetry  (0.3)  are  the  estimated  variations  due  to  ridge 
segmentation,  while  errors  in  AT  are  ±25°C,  as  defined  according  model  A  results. 


observed  MBA  and  bathymetry  along  isochrons  are  due  primarily 
to  the  crustal  structure  that  was  frozen  into  the  lithosphere  at  the 
time  of  accretion. 

Our  linear  melting  function  yields  an  essentially  linear  relation 
between  AT  and  MBA  and  bathymetry.  This  relation  is  derived 
empirically  by  a  least  squares  fit  between  theoretical  values  of 
MBA  and  bathymetry  and  corresponding  values  of  base  layer  AT 
for  model  A  calculations.  Only  points  further  than  ~80  km  from 
transform  offsets  are  used  in  the  fit.  The  dependence  of  AT  on 
MBA  is  found  to  be 

AT=-0.576AMBA  (4) 


with  a  <  3°C  standard  deviation  misfit  to  model  calculations.  The 
dependence  of  AT  on  depth  anomaly  A//  is  found  to  be 

AT=48.3AT/  (5) 

with  a  <  5°C  standard  deviation  misfit  to  model  calculations. 

Using  the  peak  mantle  Bouguer  and  bathymetry  anomalies 
along  each  isochron,  we  derive  peak  temperature  anomalies 
beneath  paleospreading  centers  (Figure  7).  Along  the  7.7  Ma 
isochron  the  observed  MBA  is  -150  mGal,  and  bathymetry 


a) 


Longitude 


Figure  8.  (a)  Peak  base  layer  AT  calculated  from  MBA  (circles)  and  bathymetry  anomalies  (triangles)  along 
selected  isochrons  are  plotted  against  the  isochron  ages.  The  solid  line  marks  AT  averaged  between  the  gravity  and 
bathymetry  calculations.  The  uncertainty  of  d:25°C  in  AT  is  the  uncertainty  estimated  from  results  for  the  present- 
day  ridge  axis,  (b)  Map  showing  the  Galapagos  hot  spot  and  the  locations  of  peak  temperature  anomalies  along  the 
present  and  paleoaxes  of  the  Gal^ipagos  spreading  center.  The  3.0-,  2.0-,  and  1.0-km  depth  intervals  are  contoured, 
and  the  ridge  axis  is  marked  as  a  bold  solid  line. 
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anomaly  is  -1.7  km;  these  anomalies  yield  a  past  temperature 
anomaly  of  ~86±25®C,  70%  greater  than  the  anomaly  along  the 
present-day  ridge  axis.  As  shown  in  Figure  8a  both  MBA 
(circles)  and  bathymetry  (triangles)  relationships  produce 
consistent  temperature  anomalies. 

Also  illustrated  in  Figure  8a  is  the  decrease  in  amplitudes  of 
the  MBA  and  bathymetry  anomalies  with  decreasing  isochron 
age.  This  behavior  indicates  that  the  peak  temperature  anomaly 
beneath  the  Galapagos  spreading  center  has  steadily  decreased 
since  7.7  Ma,  when  the  hot  spot  was  at  or  near  the  ridge  axis.  As 
the  hot  spot  migrated  southwest  away  from  the  ridge  axis  beneath 
the  Nazca  Plate,  the  maximum  in  AT"  decreased  and  the  axial 
position  of  the  peak  in  AT  moved  westward  at  approximately  the 
same  rate  as  the  westward  velocity  component  of  the  hot  spot 
with  respect  to  the  Cocos  Plate  (Figure  8b).  If  we  assume  that  the 
temperature  anomaly  of  the  hot  spot  has  remained  constant  over 
the  past  -8  m.y.,  the  above  results  provide  evidence  that  the 
amplitude  of  temperature  anomaly  beneath  the  Galapagos  ridge 
axis  is  a  function  of  the  distance  separating  the  hot  spot  and  ridge 
axis.  Such  a  dependence  may  reflect  the  cooling  of  plume 
material  as  it  migrates  from  the  Galapagos  hot  spot  to  the  ridge 
axis  and  may  provide  importance  constraints  on  the  mechanisms 
of  heat  transfer  between  hot  spots  and  nearby  ridges. 

Condusions 

The  2-D  pattern  of  the  mantle  Bouguer  and  bathymetry 
anomalies  reflect  temperature-dependent  density  structure 
imposed  by  the  Galapagos  hot  spot.  Correlation  of  MBA  and 
bathymetry  with  geochemical  anomalies  supports  the  notion  of 
mixing  of  a  hot,  enriched  plume  with  the  cooler,  depleted  upper 
mantle.  Profiles  of  mantle  Bouguer  gravity  anomalies  taken 
along  isochrons  of  ages  0.0-7.67  Ma  indicate  that  long- 
wavelength  topography  is  isostatically  compensated  by  density 
structure  in  the  crust  and  upper  100  km  of  the  mantle.  To  account 
for  the  -90  mGal  along-axis  decrease  in  MBA  and  the  -1.1  km 
decrease  in  depth,  our  models  require  a  subaxial  temperature 
anomaly  of  50±25°C  and  an  associated  crustal  thickness  increase 
of  3±1  km.  Mantle  temperatures  decrease  dramatically  from  the 
hot  spot  to  the  ridge  axis  but  decrease  much  more  gradually  along 
axis  with  a  lateral  temperature  gradient  10-20  times  less.  This 
contrast  places  important  constraints  on  hot  spot-to-ridge  and 
along-ridge  heat  transport. 

From  the  crustal  isochron  of  age  7.7  M a  to  the  present-day 
axis,  the  along-isochron  amplitudes  of  MBA  decrease  from  -150 
to  -90  mGal.  The  corresponding  along-isochron  bathymetry 
anomalies  decrease  from  -1.7  to  -1.1  km.  These  MBA  and 
bathymetry  anomalies  indicate  that  the  axial  temperature  anomaly 
was  70%  hotter  in  the  past  (86±25®C)  and  has  steadily  decreased 
to  50±25°C  as  the  ridge  axis  migrated  away  from  the  Galapagos 
hot  spot.  The  simplest  explanation  for  this  apparent  decrease  in 
the  mantle  anomaly  beneath  the  Galapagos  spreading  center  since 
7.7  Ma  is  that  the  ridge  axis  temperature  structure  depends  on  the 
distance  separating  the  hot  spot  and  ridge  axis.  These 
conclusions  point  to  the  need  for  further  experimental  and 
numerical  investigations  to  better  understand  the  dynamic 
interaction  between  the  Galapagos  spreading  center  and  hot  spot 
and  the  effects  of  such  interactive  processes  on  the  internal 
structure  of  the  oceanic  lithosphere. 
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ABSTRACT 

We  analyzed  bathymetric  and  gravity  anomalies  along  present 
and  paieoaxes  of  oceanic  spreading  centers  influenced  by  the  Ice¬ 
land,  Azores,  Galapagos,  Tristan,  and  Easter  hotspots.  Residual 
bathymetry  (up  to  4.7  km)  and  mantle  Bouguer  gravity  (up  to  —340 
km)  anomalies  are  maximum  at  on-axis  hotspots  and  decrease  with 
increasing  ridge-hotspot  separation  distance  (D),  until  becoming 
insignificant  at  D  —500  km.  Along-isochron  widths  of  bathymetric 
anomalies  (up  to  2700  km)  depend  inversely  on  paleo-spreading 
rate,  reflecting  the  extent  to  which  plume  material  will  flow  along 
axis  before  being  swept  away  by  the  spreading  lithosphere.  Flux 
balance  arguments  suggest  that  the  five  hotspots  feed  material  to 
ridges  with  comparable  fluxes  of  — 2J2  x  10*  km^/m.y.  Assuming  that 
the  amplitudes  of  these  geophysical  anomalies  reflect  temperature-de- 
pendent  crustal  thickness  and  mantle  density  variations,  we  suggest 
that  ridge  temperature  anomalies  are  maximum  (150-225  °C)  when 
plumes  are  ridge  centered  and  decrease  with  increasing  ridge-hotspot 
distance  due  to  cooling  of  the  ridgeward-migrating  plume  material. 

INTRODUCTION 

When  mantle  plumes  rise  near  oceanic  spreading  centers,  they 
generate  not  only  near-ridge  hotspots,  but  also  melt  anomalies  at 
the  axis  of  the  nearby  ridges  (e.g.,  Morgan,  1978).  Direct  evidence 
that  near-ridge  plumes  divert  toward  and  feed  ridges  is  the  ocean- 


island  basalt  (OIB)  geochemical  signature  in  ridge  basalts  (e.g..  Hart 
et  al.,  1973).  Furthermore,  along-axis  gradients  in  the  strength  of  OIB 
signatures  and  in  topography  (e.g.,  Vogt,  1976;  Schilling,  1991)  indicate 
that  once  a  plume  reaches  a  ridge,  it  spreads  laterally  along  axis. 

Previous  studies  of  ridge-plume  interactions  have  focused  pri¬ 
marily  on  present-day  spreading  centers.  Ito  and  Lin  (1995),  how¬ 
ever,  demonstrated  that  70%-75%  of  off-axis  bathymetric  and  grav¬ 
ity  anomalies  of  the  Cocos  plate  can  be  attributed  to  the  anomalous 
crustal  thicknesses  generated  at  the  paleo-Galapagos  ridge  axis.  We 
attributed  long-wavelength  (>200  km)  variations  in  bathymetry  and 
gravity  along  crustal  isochrons  to  temperature  conditions  beneath 
the  hotspot-influenced  ridge  axis  at  the  time  the  crust  was  created. 

In  this  study  we  investigated  the  evolution  of  five  prominent 
plume-ridge  systems  over  wide  ranges  in  ridge-hotspot  separation 
distance  and  spreading  rate.  The  results  of  this  study  provide  ob¬ 
servational  constraints  on  the  amplitudes  of  along-isochron  bathy¬ 
metric  and  gravity  anomalies  as  they  depend  on  ridge-hotspot  sep¬ 
aration  distance,  and  along-isochron  widths  of  bathymetric  anomalies 
as  they  depend  on  ridge  spreading  rate. 

ALONG-ISOCHRON  BATHYMETRIC  AND  GRAVITY 
ANOMALIES 

Iceland,  Azores,  Tristan,  Galapagos,  and  Easter  (Fig.  1)  are  the 
five  hotspots  that  impose  the  most  prominent  bathymetric  and  geo- 


Figure  1.  Regional  bathymetric  maps  (merca- 
tor  projections)  of  five  prominent  hotspot- 
ridge  systems:  Iceland,  Azores,  and  Tristan, 
near  Mid-Atlantic  Ridge;  Galapagos,  near  Ga¬ 
lapagos  spreading  center;  and  Easter,  near 
East  Pacific  Rise.  EtopoS  (Earth  topography  at 
5  minute  grid  spacing.  National  Geophysical 
Data  Center  Report  MGG-5)  bathymetry 
points  within  5  min  of  ship  data  points  were 
omitted  before  gridding  at  5  min  grid  spacing. 
Circles  mark  present-day  locations  of  hot¬ 
spots;  solid  lines  mark  ridge  axes  and  off-axis 
isochrons  along  which  data  profiles  were 
taken.  To  exclude  sea  floor  affected  by  off- 
axis  volcanism  we  used  isochrons  of  ages 
0-^  Ma  for  Iceland  and  0-25  Ma  for  Azores 
on  North  American  plate;  0-8  Ma  for  Galapa¬ 
gos  on  Cocos  plate;  and  0-20  Ma  for  Easter  on 
Pacific  plate.  For  Tristan,  we  used  isochrons 
of  ages  0-70  Ma  on  South  American  plate  and 
ages  80-110  Ma  on  African  plate  because  hot¬ 
spot  crossed  from  South  American  to  African 
plate  at  —80  Ma  (O’Connor  and  Duncan,  1990). 
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chemical  anomalies  observed  at  nearby  oceanic  spreading  centers 
(Hart  et  aL,  1973;  Hamelin  et  al.,  1984;  Schilling,  1985).  Encom¬ 
passing  each  of  the  five  systems,  we  obtained  shipboard  bathymetric 
data  from  the  National  Geophysical  Data  Center  (NGDC)  and 
Lamont-Doherty  Earth  Observatory  (LDEO),  and  gridded  bathym¬ 
etry  from  NGDC.  To  derive  residual  bathymetry,  we  first  corrected  the 
raw  data  for  isostatic  effects  of  sediment  loading  and  then  subtracted 
predicted  depths  of  a  cooling  mantle  half  space  (Carlson  and  Johnson, 
1994).  Sediment  thicknesses  were  obtained  from  the  LDEO  database 
(A,  Cazenave,  Centre  National  d’Etudes  Spatiales,  Toulouse,  France), 
and  density  contrasts  between  the  sediments  and  mantle,  and  mantle 
and  water  were  assumed  to  be  1600  kg/m^  and  2300  kg/m^,  respectively. 

Free-air  gravity  data  were  taken  from  the  ship  surveys  and  the 
satellite  altimetry-derived  gravity  grid  of  Sandwell  and  Smith  (1992). 
To  isolate  the  effects  of  sub-sea-floor  density  structure,  we  gener¬ 
ated  mantle  Bouguer  anomalies  by  subtracting  from  the  free-air 
gravity  the  attractions  of  the  sea-floor-water  (density  contrast,  Ap  = 
1800  kg/m^)  and  crust-mantle  (Ap  =  500  kg/m^)  interfaces  using  raw 
bathymetry,  and  assuming  a  crust  of  uniform  thickness  (6.5  km) 
(e.g.,  Kuo  and  Forsyth,  1988). 

Coordinates  of  present-day  ridge  axes  and  crustal  isochrons 
were  defined  by  using  plate  boundary  and  age  data  of  Muller  et  al. 
(1993a).  Because  our  focus  was  on  anomalies  generated  at  the  axes 
of  spreading  centers,  we  considered  only  data  from  sea  floor  unaf¬ 
fected  by  off-axis  volcanism,  as  detailed  in  the  Figure  1  caption. 
From  our  residual  bathymetry  and  mantle  Bouguer  grids,  we  then 
extracted  along-isochron  profiles  (Fig.  2). 

ANOMALY  AMPLITUDES  VS.  PALEORIDGE-HOTSPOT 
DISTANCE 

To  determine  hotspot  locations  relative  to  paleo-spreading 
centers,  we  assumed  that  the  hotspots  were  stationary  with  respect 
to  each  other  and  used  plate-reconstruction  poles  (Lonsdale  1988; 


Figure  2.  (A)  Residual  bathymetry  (RB)  and  (B)  mantle  Bouguer  anomaly 
(MBA)  profiles  along  six  example  isochrons  of  Tristan  system.  Shaded 
parts  mark  long-wavelength  signals  we  attribute  to  hotspots.  W  de¬ 
fines  along-isochron  width  over  which  long-wavelength  topographic 
swells  are  shallower  than  depths  predicted  by  cooling  half-space  ref¬ 
erence  model.  AHfl  and  AMBA  are  maximum  amplitudes  along  each 
profile.  Decrease  in  amplitudes  with  decreasing  isochron  age  coincide 
with  migration  of  Tristan  hotspot  away  from  ridge  axis  since  ~80  Ma 
when  it  was  ridge  centered. 
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Miiller  et  al.,  1993b)  to  rotate  isochrons  with  respect  to  the  hotspots 
back  to  their  positions  at  the  time  of  accretion.  We  then  measured 
distances  between  the  paleo-ridge  axes  and  hotspot  centers,  which 
we  took  to  be  the  locations  of  most  recent  volcanism. 

The  along-isochron  variations  in  residual  bathymetric  (ARB) 
and  mantle  Bouguer  anomalies  {AMBA)  display  a  decrease  with 
increasing  paleoridge-hotspot  distance  (D,  Fig.  3).  The  on-ridge  hot¬ 
spot  cases  (D  <50  km)  for  the  Tristan  system  (80-90  Ma  isochrons) 
and  the  Iceland  system  (0-30  Ma  isochrons)  display  the  highest  ARB 
(3.5- 4.7  km)  and  most  negative  AMBA  (-250  to  -340  mgal),  which 
are  approximately  twice  those  of  ridge-centered  cases  for  the  Ga¬ 
lapagos  and  Azores  systems.  At  D  ~  500  km,  the  hotspot  signals 
become  very  weak  and  in  the  case  of  Tristan,  become  indistinguish¬ 
able  from  normal  ridge-segmentation-related  variations.  The  indi¬ 
vidual  Galapagos  and  Tristan  systems  show  a  decrease  in  ARB  and 
AMBA  with  increasing  D,  whereas  the  Azores  system  is  more  com¬ 
plex  and  the  Easter  trend  is  very  weak.  The  predominant  decrease 
of  ARB  and  AMBA  with  increasing  D  is  consistent  with  Schilling’s 
(1985)  study  of  present-day  ridge-axis  bathymetry. 

ANOMALY  WIDTHS  VS.  PALEO-SPREADING  RATE 

Whereas  amplitudes  of  ARB  and  AMBA  are  functions  of  ridge- 
hotspot  distance,  along-isochron  widths  (W)  of  the  bathymetric 
anomalies  (see  Fig.  2  caption)  depend  primarily  on  the  full  spread¬ 
ing  rate  (U)  at  the  time  of  crustal  accretion.  The  maximum  values 
of  W  are  found  along  the  slowest-spreading  Mid-Atlantic  Ridge 
near  Iceland  (2700  km.  Fig.  4A);  these  values  are  comparable  to  the 
along-axis  extent  of  the  helium  isotope  anomaly,  but  are  a  factor  of 
two  greater  than  the  widths  of  rare-earth-element  anomalies  (Schill¬ 
ing,  1986).  Values  of  W  decrease  with  increasing  U  to  a  minimum 
along  the  fast-spreading  East  Pacific  Rise. 

The  observed  dependence  of  W  on  spreading  rate  lends  strong 
support  to  previous  notions  of  along-axis  plume  material  flow  (Vogt, 
1976;  Schilling,  1985).  Similarly  to  Schilling  (1991),  we  estimate  that 
the  flux  of  plume  material  feeding  the  ridge  (Q)  is  eventually  carried 
away  by  the  spreading  lithospheric  plates  (Fig.  4A,  inset),  such  that 

hUW 

Q=l  P(y)hUdy  =  —,  (1) 

where  y  is  the  along-axis  coordinate,  h  is  the  thickness  of  the  fully 
developed  lithosphere  (assumed  to  be  80  km),  and  P{y)  is  the  per¬ 
centage  of  accreted  lithosphere  derived  from  the  plume  material 
assumed  to  decrease  linearly  from  1  aty  =  0  to  0  aty  =  ±W12. 

We  treat  the  hotspot  to  ridge  flow  as  a  simple  laminar  flow 
problem  in  which  the  lithospheric  drag  opposes  the  ridgeward  flow 
of  plume  material.  The  channel  connecting  the  ridge  and  hotspot 
has  a  characteristic  width  and  thickness  w-2  (see  Fig.  4,  A  and  B, 
insets).  Therefore,  the  net  flux  from  the  hotspot  to  the  ridge  is 

(2) 

where  K  is  the  average  ridgeward  velocity  of  plume  flow,  is  the 

ridgeward  flux,  and  is  the  opposing  plate-driven  flux.  (Com¬ 

bining  equations  1  and  2  yields  the  dependence  of  W  on  spreading  rate. 


The  solid  curve  in  Figure  4A  is  that  predicted  for  assumed  values  of 
K  —  70  km/m.y.  and  w^W2  =*  3  X  10^  km^,  which  yields  a  root-mean- 
square  misfit  to  the  data  of  500  km.  Similar  misfits  are  achieved  for 
V  =  30-100  km/m.y.  and  corresponding  values  of  W{W2  of  8-2  X  10"* 
km^.  These  results  suggest  that  the  ridgeward  fluxes  from  the  five 
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Ridge-Hot  Spot  Distance,  D  (km) 

Figure  3.  Along-isochron  amplitudes  of  (A)  ARS  and  (B)  AMBA  plotted 
against  distance  between  paleo-ridge  axes  and  hotspots  at  times  cor¬ 
responding  to  isochron  ages.  Asterisks  mark  present-day  ridge-axis 
anomalies;  solid  lines  are  those  that  best  fit  all  data.  Uncertainties  in 
ARB  and  AMBA  are  segmentation-scale  variations  (Lin  and  Phipps 
Morgan,  1992)  that  are  independent  of  larger  wavelength  hotspot  sig¬ 
nals.  Uncertainties  in  D  reflect  uncertainties  in  isochron  ages  and  in 
plate  motion  relative  to  hotspots  (Cande  and  Kent,  1992). 

hotspots  are  comparable,  the  average  value  being  —2.2  x  10^  km^l 
m.y.  Increasing  or  decreasing  ^v,^V2Kby  1  x  10^  km^/m.y.  increases 
the  data  misfit  by  a  factor  of  two. 

Our  theoretical  relation  between  W  and  U  is  based  on  one 
end-member  scenario  in  which  lateral  spreading  of  plume  material 
beneath  ridges  is  strictly  ridge  parallel.  A  numerical  study  that  con¬ 
siders  both  ridge-parallel  and  ridge-perpendicular  spreading  of 
plumes  beneath  ridges  (Feighner  et  al.,  1995)  may  represent  the 
other  end  member;  it  thus  predicts  that  Wis  proportional  to  (QIU)^^ 
rather  than  (Q/U),  as  does  our  model. 

PALEO-RIDGE-AXIS  TEMPERATURE  ANOMALIES 

We  show  here  how  the  amplitudes  of  and  dMBA  may  reflect 
the  temperature  imomalies  beneath  the  paleo-ridge  axes.  We  assume 
that  ARB  and  AMBA  arise  from  crustal-thickness  and  mantle-density 
anomalies,  both  of  which  depend  on  the  ridge-axis  temperature  anom¬ 
aly  at  the  time  of  crustal  accretion.  ARB  and  AMBA  can  be  related  to 
a  hoLspot-induced  mantle  temperature  anomaly  (AT)  using  the  model 
of  Ito  and  Lm  (1995),  which  considered  changes  in  mantle  density  by 
thermal  expansion,  and  in  cmstal  thickness  by  increased  decompres¬ 
sion  melting.  Assuming  passive  mantle  upwclling,  Ito  and  Lin  (1995) 
impo.sed  temperature  anomalies  below  the  melting  zone  and  then  com¬ 
bined  the  effects  of  crustal-thickness  and  mantle-density  variations  to 
yield  theoretical  isostatic  bathymetric  variations  and  AMBA. 

Applying  this  method  for  ranges  of  imposed  temperature  anom¬ 
alies  and  model  spreading  rates,  we  derive  the  empirical  relations: 

AT=(0.11C/  +  35.3)A/?5,  (4) 


Figure  4.  A:  Along-isochron  widths  of  residual  bathymetric  anomalies 
vs.  full  spreading  rates  during  times  corresponding  to  isochron  ages. 
Diamonds — Iceland,  triangles — ^Azores,  squares — Galapagos,  light 
blue  circles — Tristan,  dark  blue  circles— -Easter.  Asterisks  mark 
present-day  ridge-axis  anomalies.  Present-day  spreading  rates  are 
from  DeMets  et  al.  (1990);  paleorates  are  from  Cande  and  Kent  (1992), 
Mayes  et  al.  (1990),  Srivastava  and  Tapscott  (1986),  and  Wilson  and  Hey 
(1995).  Note  that  Iceland,  Azores,  Galapagos,  and  Easter  plot  in  tight 
groupings,  each  defining  narrow  range  in  W  (-^600  km)  and  U  (:^20 
km/m.y.).  Tristan  anomalies  encompass  wider  range  in  W,  reflecting 
secondary  dependence  on  D  (best  fitting  line  is  IV  =  1690  —  3.3D; 
root-mean-square  misfit  is  250  km).  Solid  curve  is  relation  derived  in 
text.  Inset  illustrates  map  view  of  plume-ridge  flow  pattern;  dot  pattern 
marks  plume  material.  B:  Along-isochron  temperature  anomalies  (de¬ 
rived  from  method  of  Ito  and  Lin,  1995)  plotted  against  D.  Solid  curve 
is  predicted  by  conductive  cooling.  Inset  illustrates  depth  cross  sec¬ 
tion  of  plume  conduit  between  hotspot  and  ridge  axis. 

and 

at  =  -{0.001717  0.45)AMSyl.  (5) 

The  dependence  of  AT  on  U  reflecLs  a  subtle  dependence  of  crustal 
thickness  on  spreading  rate  that  is  consistent  with  calculations  of  Su 
et  al.  (1994).  For  AT  =  100  X  and  U  =  20-100  km/m.y.,  for  example, 
we  predict  correxSponding  values  of  crustal  thlckenmg  of  9-4.5  km. 

Temperature  anomalies  derived  accordingly  from  the  ob.scrved 
ARB  and  AMBA  are  maximum  for  the  on-ridge  cases  (150-225  "C), 
and  decrease  to  near  zero  for/)  --  500  km  (Fig.  4B).  Such  a  behavior 
can  be  interpreted  as  the  cooling  trend  of  plume  material  as  it 
migrates  from  hotspot  centers  to  nearby  ridges,  the  ridge-centered 
cases  reflecting  the  temperature  anomaly  of  the  hotspot  itself. 

As  plume  material  migrates  from  a  hotspot  center  to  a  ridge,  it 
conducts  heat  to  the  surrounding  mantle  (see  Fig.  4,  A  and  B,  in¬ 
sets).  Assuming  that  the  amount  of  heat  conducted  in  the  direction 
of  plume  flow  is  negligible,  the  heat  balance  equation  is 
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where  T  is  the  average  temperature  of  the  plume  conduit,  p  and 
are  the  density  and  heat  capacity  of  the  plume  material,  respectively, 
and  qy  and  are  the  components  of  conductive  heat  flow  out  of  the 
conduit  walls.  If  we  assume  that  heat  loss  occurs  through  a  thermal 
boundary  layer  surrounding  the  plume  channel  with  thickness  8, 


k^T 


(7) 


where  k  is  the  mantle  thermal  conductivity  (3  W^m"^  and  8 
is  defined  such  that  q  =  30-100  mW/m^  comparable  to  heat-flux 
values  on  the  sea  floor.  If  it  is  assumed  that  the  gradients  dqyidy  and 
dqjdz  are  proportional  to  1/iVj  and  lAv2,  respectively,  and  Q  ~ 
^1^2^  ^  ^  W),  the  combination  of  equations  6  and  7  yields 


BT  2k 


(8) 


where  k  =  kfpCp  is  thermal  diflhisivity  (10  ®  m^/s).  Integrating  with 
respect  to  x  from  0  to  D  yields 


AT  =  AT^exp 


2k 

(w^i  +  w^D 


(9) 


where  AT^  is  the  temperature  anomaly  at  the  hotspot  center. 

Taking  AT^  =  100  ®C,  (wj  +  Wj)  =  400  km,  and  0  =  2.2  x  10^ 
km^/m.y.  as  consistent  with  the  observed  W  vs.  U  trend  above,  we 
produce  a  theoretical  curve  (Fig.  4B)  that  effectively  matches  the 
inferred  temperature  anomalies  for  Z)  >  50  km.  For  D  ^  100  km, 
the  Iceland  and  Tristan  points  lie  significantly  higher  than  the  the¬ 
oretical  curve.  This  mismatch  may  be  because  (1)  the  Iceland  and 
early  Tristan  plumes  are  hotter  than  the  other  hotspots  and/or  (2) 
latent  heat  loss  due  to  melting  at  the  hotspot  centers  rapidly  cools 
the  plume  before  it  migrates  to  nearby  ridges  in  the  off-axis  cases. 
For  D  500  km,  AT  is  small  enough  that  its  effects  on  bJRB  and 
bMBA  are  negligible,  even  though  the  plume  may  still  be  feeding  the 
ridge.  Consequently,  the  geochemical  signal  may  persist  to  a  ridge- 
hotspot  distance  of  up  to  850  km  (Schilling  et  al.,  1985),  long  after 
the  signals  in  bKB  and  AA/B^  have  disappeared. 


CONCLUSIONS 

Along-isochron  variations  in  residual  bathymetry  and  mantle 
Bouguer  gravity  reflect  the  influence  of  hotspots  on  paleoaxes  of 
nearby  spreading  centers.  The  amplitudes  of  dong-isochron  anom¬ 
alies  for  the  five  prominent  plume-ridge  systems  reach  a  maximum 
of  4.7  km  for  bKB  and  —340  mgal  for  ^MBA  and  decrease  with 
increasing  paleoridge-hotspot  distance.  The  along-isochron  widths 
(0-2700  km),  however,  depend  inversely  on  paleo-spreading  rate. 
Whereas  the  widths  of  ABB  reflect  the  balance  between  ridgeward 
plume  flux  and  lithospheric  accretion,  the  amplitudes  of  ABB  and 
MBA  reflect  paleoaxial  temperature  anomalies  that  decrease  as 
the  plume  material  cools  along  its  lateral  migration  to  nearby  ridges. 
The  five  hotspots  appear  to  deliver  material  to  ridges  with  compa¬ 
rable  fluxes  of  ^2.2  X  10^  km^/m.y,  and  produce  excess  mantle 
temperature  anomalies  of  50  to  225  ®C  that  influence  ridge-axis 
structure  to  a  maximum  ridge-hotspot  distance  of  --500  km. 
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Abstract 

We  investigate  the  dynamics  of  mantle  flow  and  melting  of  a  ridge-centered  plume  with 
three-dimensional  variable-viscosity  numerical  models,  focusing  on  three  buoyancy 
sources:  temperature,  melt  depletion,  and  melt  retention.  The  width  W,  to  which  a  plume 
spreads  along  a  ridge  axis,  depends  on  plume  volume  flux  Q,  full  spreading  rate  U, 
buoyancy  number  B,  and  ambient/plume  viscosity  contrast  y.  When  all  melting  effects  are 
considered,  our  numerical  results  are  best  parameterized  by  'W=2.2)l{Q/U)^^^{B'f)^-^^. 
Thermal  buoyancy  is  first  order  in  controlling  along-axis  plume  spreading  while  latent  heat 
loss  due  to  melting,  and  depletion  and  retention  buoyancy  forces  contribute  second-order 
effects.  We  propose  two  end-member  models  for  the  Iceland  plume  beneath  the  Mid- 
Atlantic  Ridge  (MAR).  The  first  has  a  broad  plume  source  with  temperature  anomaly  ATp 
of  75°C,  radius  a  of  300  km,  and  Q  of  1 .2x10^  km^/my.  The  second  is  of  a  narrower  and 
hotter  plume  source  with  ATp  of  170°C,  radius  of  60  km,  and  Q  of  2.1x10^  km^/my.  The 
broad  plume  source  predicts  successfully  the  observed  seismic  crustal  thickness, 
topographic,  and  gravity  anomalies  along  the  MAR,  but  predicts  an  along-axis  geochemical 
plume  width  substantially  broader  than  that  suggested  by  the  observed  Sr^^/Sr^®  anomaly. 
The  narrow  plume  source  model  predicts  successfully  the  total  excess  crustal  production 
rate  along  the  MAR  (2.5xl0^km3/my)  and  a  geochemical  width  consistent  with  that  of  the 
Sr^^/Sr^^  anomaly,  but  it  requires  substantial  along-axis  melt  transport  to  explain  the 
observed  along-axis  variations  in  crustal  thickness,  bathymetry,  and  gravity.  Calculations 
suggest  that  lateral  plume  dispersion  may  be  radially  symmetric  rather  than  channeled  along 
the  ridge  axis  and  that  the  topographic  swell,  which  is  elongated  along  the  Reykjanes 
Ridge,  may  be  due  to  rapid  off-axis  subsidence  associated  with  lithospheric  cooling 
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superimposed  on  a  broader  hotspot  swell.  The  two  plume  source  models  predict  seismic 
P-wave  velocity  reductions  of  0.5-2%  in  the  center  of  the  plume,  producing  travel  time 
delays  of  0.2-1 .2  s.  Predicted  P-wave  delay-times  for  the  narrow  plume  source  model  are 
more  consistent  with  recent  seismic  observations  beneath  Iceland,  suggesting  that  this 
model  may  be  more  representative  of  the  Iceland  plume. 


1.  Introduction 

Centered  on  the  Mid-Atlantic  Ridge  (MAR),  the  Iceland  hotspot  is  the  largest  melt 
anomaly  throughout  the  world’s  mid-ocean  ridge  system  and  is  among  the  large  oceanic 
igneous  provinces  [1],  The  idea  that  Iceland  marks  a  mantle  convection  plume  rising 
beneath  the  MAR  has  become  well  established  since  its  original  conception  in  the  early 
1970’s  (e.g.  ref.  [2,  3,  4]).  The  broad  topographic  swell  (Fig.  1)  and  correlated  along- 
spreading-axis  geochemical  anomalies  indicate  that  the  plume  rises  beneath  Iceland  and 
spreads  laterally  along  the  ridge  axis  [4,  5].  Such  along-axis  spreading  of  a  mantle  plume 
feeding  a  ridge  axis  may  also  explain  topographic  and  geochemical  anomalies  affected  by 
other  near  ridge-axis  hotspots  (e.g.  ref  [6,  7,  8]) — many  of  which  may  have  contributed 
substantially  to  the  earth’s  heat  and  magmatic  budget  throughout  geologic  history. 

While  the  original  concept  that  plumes  feed  and  spread  along  nearby  ridges  was 
proposed  two  decades  ago,  only  recently  have  the  fluid  dynamic  aspects  been  investigated 
quantitatively.  Recent  numerical  and  laboratory  tank  experiments  have  shown  that  the 
width  W,  over  which  a  plume  spreads  along  axis,  increases  with  plume  volume  flux  Q,  and 
decreases  with  plate  full-spreading  rate  U  [9,  10,  11].  Such  studies  are  important  in 
revealing  the  pertinent  physical  processes  governing  plume-ridge  interactions  and  in  placing 
theoretical  constraints  on  properties  of  mantle  plumes  such  as  temperature  anomaly,  size, 
and  volume  flux. 

Two  potentially  important  sources  of  buoyancy,  however,  have  not  been  considered  in 
previous  plume-ridge  studies.  These  are  melt  depletion,  which  lowers  the  Fe/Mg  ratio  in 
the  residual  mantle  and  thus  reduces  its  density  [12],  and  melt  retention  in  the  mantle, 
which  also  reduces  mantle  bulk  density  (e.g.  [13,  14,  15]).  It  has  been  proposed  that  melt 
depletion  may  be  primary  in  driving  spreading  of  intraplate  plumes  beneath  the  lithosphere 
[16].  It  has  also  been  proposed  that  both  melt  retention  buoyancy  and  depletion  buoyancy 
may  contribute  significantly  to  along-axis  variations  in  mantle  flow  and  crustal  thickness 
beneath  normal  mid-ocean  ridges  [17,  18]. 
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The  objectives  for  this  study  are  two  fold.  First  we  investigate  numerically  the  effects 
of  thermal-  and  melting-related  buoyancy  forces  on  along-axis  spreading  of  ridge-centered 
plumes.  We  use  three-dimensional  (3D),  variable  viscosity,  numerical  models  to  simulate  a 
buoyant  plume  rising  beneath  spreading  plates  and  systematically  test  the  effects  of  thermal, 
melt  depletion,  and  melt  retention  buoyancy  forces.  Our  second  objective  is  to  constrain 
the  temperature  anomaly,  dimension,  and  volume  flux  of  the  Iceland  plume  by  comparing 
theoretical  predictions  with  observed  variations  in  seismic  crustal  thickness,  topography, 
gravity,  and  geochemistry  on  Iceland  and  along  the  Mid-Atlantic  Ridge.  We  propose  two 
end-member  models  for  the  mantle  plume  source  beneath  Iceland  to  explain  the 
observations,  and  discuss  their  implications  on  basalt  geochemistry,  melt  migration,  and 
seismic  velocity  variations  along  the  Mid-Atlantic  Ridge  axis. 

2.  Governing  equations 

To  model  mantle  flow  of  a  plume-ridge  system  in  the  upper-mantle,  we  treat  the  mantle 
as  a  fluid  of  zero  Reynolds  number  and  infinite  Prandtl  number.  The  3D  stress  tensor,  T,  is 

defined  according  to 

x=2t]{Tr,p)  e -pi  (1) 

where  I  is  the  identity  matrix  and  T]  is  viscosity  which  depends  on  real  temperature  Tr  and 
hydrostatic  pressure  p.  The  strain  rate  tensor  e  depends  on  spatial  derivatives  of  mantle 
flow  rate  u  according  to  e  =  U2(uij  +  uji  ).  The  equilibrium  equations  include 
conservation  of  mass 

V*«  =  0,  (2) 

momentum 

V.T  =  -Ap(r,X,(^>)g£,  (3) 

and  energy 

dT  xi'T  rA\ 

—  =k\T-u*\T - M  (4) 

ot  Cp 

(see  Table  1  for  definition  of  variables).  Eq.  (2)  satisfies  the  Boussinesq  approximation 
and  neglects  dilational  flow  due  to  the  extraction  of  melt  which  is  likely  to  be  small  [19]. 
Eq.  (3)  balances  viscous  stresses  with  the  body  force  due  to  density  variations  which 
depend  on  potential  temperature  T,  melt  depletion  X,  and  mantle  porosity  0,  according  to, 
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(5) 


Ap  =  -Po 


aT  +  (3X  + 


Po  Pm  0 
Po  ) 


Eq.  (4)  balances  energy  transfer  associated  with  heat  conduction,  heat  advection,  and  latent 
heat  loss  due  to  melting.  Melt  depletion  is  governed  by 


dt 


(6) 


where  M  is  melt  fraction  and  M  = 

dt 

To  estimate  the  distribution  of  porosity  (p,  we  assume  that  melt  migrates  vertically 
through  the  mantle  at  a  melt-mantle  velocity  contrast  {(O-w)  as  governed  by  Darcy's  flow 
law. 


(pio)  -  w)  = 


{Po  Pm 


Pn 


(7) 


^2^2 

Permeability  K  depends  on  grain  size  b  according  to  AT  =  — Finally,  the  rate  of  melt 

12k 

percolation  is  assumed  to  be  equivalent  to  the  rate  at  which  melt  is  generated  such  that 


(p{z)a){z)  =  [ Mdz. 

Pm  ^ 


(8) 


3.  Numerical  method  and  boundary  conditions 

To  solve  the  above  equations,  we  use  a  Cartesian  numerical  code  presented  by  Gable 
[20,  21],  Time  integration  is  achieved  by  iterating  through  discrete  time  steps,  during  each 
of  which  we  solve  for  mantle  flow,  mantle  potential  temperature,  and  melt  depletion.  In 
solving  the  dimensionless  forms  of  the  flow  equations  (Eqs.  1-4),  horizontal  derivatives 
are  expressed  in  terms  of  their  Fourier  components  while  vertical  derivatives  are  expressed 
as  finite  difference  approximations.  We  then  invert  for  horizontal  and  vertical  components 
of  velocities  and  stresses  using  a  standard  relaxation  method. 

The  dimensionless  form  of  Eq.  (3)  is 
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(9) 


3  /  \ 

V.f  =  aT^T+pX  +  ^-^~^'^  <j)\, 

V  Po  ) 

where  primes  denote  dimensionless  variables.  The  body  force  (right  hand  side  of  Eq.  (9)) 

is  the  sum  of  three  terms:  the  first  term,  which  scales  with  T,  is  a  Rayleigh  number. 


ocTfj\  the  second  term,  which  scales  with  X,  is  a  melt 


depletion  Rayleigh 


number,  P; 


and  the  third  term,  which  scales  with  0,  is  a  melt  retention 


Rayleigh  number,  Ra(^  PoS^  Po — Pni  Assumed  values  for  /?  and  — — are 

^Po  \  Po  )  \  Po  j 

0.06  [12,  16]  and  0.121  [14,  18],  respectively.  Consequently,  depleting  the  mantle  by 
25%  yields  a  density  reduction  equivalent  to  heating  the  mantle  by  440°C,  while  a  melt 
porosity  of  3%  yields  a  density  reduction  equivalent  to  heating  the  mantle  by  107°C. 

We  assume  that  mantle  viscosity  varies  with  real  temperature  Tr  and  pressure 
according  to 


P=Po  exp 


E  +  pV 
RTr 


E  +  p^g(0.5D)V] 

RTro  J 


(10) 


where  reference  viscosity  r\o  is  defined  as  the  mantle  viscosity  for  T=To  and  z=0.5D;  Tr 
in  Kelvin  is  {T  +  0.6z  +  273),  where  the  term  0.6z  takes  into  account  the  adiabatic  gradient; 
and  Tro  is  the  real  temperature  value  of  Tq.  To  approximate  numerically  the  effects  of  non- 
Newtonian  rheology,  we  use  reduced  values  of  activation  energy  E  and  activation  volume 
V  [22]  (Table  1).  Because  lateral  variations  in  viscosity  introduce  nonlinearity  to  the  above 
flow  equations,  we  linearized  the  equations  by  introducing  additional  body  force  terms  [20, 
23].  The  nonlinear  terms  and  solutions  were  then  updated  upon  successive  iterations  until 
solutions  converged  to  our  specified  limit.  We  found  that  a  convergence  criterion  of  0. 1- 
0.5%  yielded  time-integrated  solutions  with  errors  of  <0.5%  while  minimizing  computing 
time.  This  computational  method  was  tested  in  2D  with  independent  finite  element 
solutions,  while  in  3D,  it  produced  solutions  within  2.6%  of  the  best-estimated 
extrapolated  solutions  of  a  benchmark  problem  of  ref.  [24]  . 
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The  final  velocity  field  is  then  used  in  the  advection  term  in  Eq.  (4)  to  solve  for  a  new 
temperature  field.  Our  energy  solver  uses  finite  differences  with  a  tensor  diffusion  scheme 
to  reduce  numerical  diffusion  which  is  intrinsic  to  finite  difference  methods  [20,  21],  The 
same  tensor  diffusion  method  is  used  to  solve  Eq.  (6)  for  the  depletion  field.  Vertical  flow 
determines  the  rate  of  decompression  melting,  comprising  the  source  terms  in  Eqs.  (4)  and 
(6).  The  melting-rate  term  in  Eq.  (4)  is  latent  heat  loss,  which  inhibits  buoyant  mantle  flow 
by  increasing  both  mantle  density  and  viscosity,  while  the  melting-rate  term  in  Eq.  (6) 
generates  low-density  depleted  mantle  residuum.  To  calculate  melting  rate  M,  we 
incorporate  the  solidus  and  liquidus  functions  of  McKenzie  and  Bickle  [25],  as  well  as  their 
functional  dependence  of  M  on  homologous  temperature  for  adiabatic  batch  melting. 

The  rate  of  melting  also  determines  the  volume  fraction  of  melt  retained  in  the  mantle  0, 
which  is  the  source  of  retention  buoyancy.  To  compute  porosity  we  combine  Eqs.  (7)  and 
(8)  and  solve  the  integral  in  Eq.  (8)  numerically  similar  to  ref.  [18].  The  grain-size- 
dependent  melt  permeability  that  we  incorporate  results  in  maximum  porosities  of  1-3% 
which  is  slightly  higher  than  the  0.1-1%  porosity  range  inferred  from  238u-230q'1^.226Ra 
disequilibria  in  Hawaiian  lavas  [26]. 

The  numerical  model  setup  is  illustrated  in  Fig.  2.  A  ridge  axis  (jc=0)  is  simulated  by 
defining  reflecting  temperature  (i.e.,  zero  heat  flux)  and  flow  (i.e.,  zero  shear  stress) 
boundary  conditions  at  the  vertical  sides,  and  setting  the  top  boundary  (z=0)  to  move  at  a 
constant  half-spreading  rate  0.5 f/.  Temperature  at  the  surface  (z=0)  is  maintained  at  0°C 
which  cools  and  thickens  a  high-viscosity  lithosphere  approximately  with  the  square  root  of 
X.  A  plume  is  introduced  by  imposing  a  columnar-shaped  temperature  anomaly  in  the 
lower  portion  of  the  box,  centered  beneath  the  ridge  axis.  The  plume  is  hottest 
(T=To+ATp)  at  its  center  and  cools  as  a  Gaussian  function  of  radial  distance  to  Tg  at  its 
radius  a.  We  exploit  the  symmetry  in  x  and  y  by  centering  the  plume  column  at  x=y=0, 
which  allows  a  quarter  plume  in  solution  space  to  represent  a  fully  circular  plume  in  virtual 
space.  In  the  lower  portion  of  the  box  (z  >  0.6D),  we  impose  the  potential  temperature  to 
be  To  everywhere  except  inside  the  plume  source.  Thus,  the  energy  equation  is  solved 
only  in  the  upper  portion  of  the  box  (0.6D  >  z  >  0). 

To  ensure  numerical  accuracy  in  the  flow  solutions,  we  set  a  non-dimensional  viscosity 
(77/770)  upper  limit  of  200  and  set  a  lower  limit  of  0.1.  The  upper  viscosity  limit  is 
sufficient  to  accurately  simulate  a  rigid  lithosphere  (i.e.,  u=U  and  v=w=0  in  the 
lithosphere),  while  the  lower  limit  allows  us  to  incorporate  the  full  viscosity  reduction  in  a 
plume  with  temperature  anomaly  of  200°C.  The  depth  dependence  of  viscosity  yields  a 
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factor  of  ~4  viscosity  increase  between  top  and  bottom  of  the  box  for  a  constant  mantle 
temperature. 

4.  Steady-state  along-axis  width  of  a  mantle  plume  head 

We  seek  here  to  quantify  the  effects  of  melting  on  mantle  flow  and  thus  the  dependence 
of  along-axis  plume  width  W  on  plume  flux  Q  and  plate  spreading  rate  U.  We  began 
numerical  experiments  with  the  steady-state  temperature  solution  of  a  ridge  without  the 
plume.  Then  after  activating  the  plume,  we  integrated  through  time  until  both  along-axis 
plume  width  and  plume  flux  converged  to  steady-state  values.  We  ran  four  sets  of 
experiments:  experimental  set  A  (Table  2a)  includes  only  thermal  buoyancy  and  omits  all 
melting  effects;  set  B  (Table  2b)  considers  only  thermal  buoyancy  but  includes  latent  heat 
loss;  set  C  (Table  2c)  includes  additional  buoyancy  from  melt  depletion;  and  set  D  (Table 
2d)  includes  additional  buoyancy  from  melt  retention. 

Fig.  2  shows  an  example  steady-state  velocity  and  temperature  field  for  a  calculation  in 
set  A  with  a  plume  source  temperature  anomaly  of  200°C  (model  5a).  Velocity  vectors 
illustrate  the  plume  rising  from  the  conduit  source  and  then  spreading  both  perpendicular  to 
and  along  the  ridge  axis  after  it  impinges  on  the  base  of  the  lithosphere.  Combined  effects 
of  thermal  buoyancy  and  reduced  plume  viscosity  result  in  a  maximum  plume  upwelling 
rate  of  244  km/my,  which  is  >20  times  that  of  the  half  spreading  rate  of  10  km/my.  The 
corresponding  average  upwelling  rate  in  the  melting  zone  (z  <110  km  )  is  85  km/my. 

Fig.  3a  shows  the  steady-state  velocity  and  mantle  density  fields  for  the  same  plume 
source  temperature  anomaly  but  with  the  additional  effects  of  latent  heat  loss  (model  5b). 
In  the  melting  region  of  the  plume  center,  potential  temperatures  are  - 1 30°C  cooler  and 
consequently  the  plume  is  65%  less  buoyant  and  3  times  more  viscous  than  the  calculation 
without  latent  heat  loss  (Fig.  2).  The  resulting  average  upwelling  rate  in  the  melting  zone  is 
50  km/my,  only  -60%  of  the  predicted  average  upwelling  rate  of  the  model  without  latent 
heat  loss  (model  5a). 

The  addition  of  melt  depletion  buoyancy  in  model  5c  generates  an  additional  -1%  lateral 
density  contrast  between  the  plume  center  and  the  mantle  beneath  normal  ridge  sections  far 
from  the  plume  (Fig.  3b).  The  resulting  average  melting-zone  upwelling  rate  is  67  km/my. 
As  material  rises  more  rapidly  in  the  plume  center,  it  spreads  more  rapidly  along  the  base  of 
the  rigid  lithosphere.  This  in  turn  inhibits  upwelling  at  radial  distances  of  100-150  km 
shown  as  negative  velocity  differences  in  Fig  3b. 
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Finally,  model  5d  considers  the  additional  buoyancy  from  melt  retention  (Fig.  3c).  The 
high  melting  rate  in  the  plume  center  results  in  a  maximum  porosity  of  2.5%,  to  reduce 
bulk  density  in  the  plume  center  by  an  additional  0.3%.  This  added  retention  buoyancy 
further  enhances  the  average  upwelling  rate  in  the  melting  zone  to  77  km/my,  which  is 
-90%  of  that  predicted  by  the  model  that  neglects  all  melting  effects  (model  5a).  Thus,  the 
added  melting-related  buoyancy  forces  approximately  balance  the  upwelling-inhibiting 
effects  of  latent  heat  loss. 

In  all  models  examined  we  find,  as  did  Ribe  et  al.  [11],  that  the  thickening  lithosphere 
does  not  channel  the  plume  preferentially  along  the  ridge  axis.  On  the  contrary,  the 
spreading  lithosphere  enhances  ridge-perpendicular  flow  by  pulling  plume  material  away 
from  the  ridge-axis,  and  actually  impedes  along-axis  flow  by  viscous  shear.  These  effects 
however  are  small — the  total  along-axis  flux  at  y=70  km  is  within  a  few  percent  of  the  total 
ridge-perpendicular  flux  at  x=70  km.  Thus,  the  rate  of  spreading  away  from  the  plume 
center  is  approximately  equal  in  all  radial  directions. 

To  determine  how  W  depends  on  Q  and  U  for  each  experimental  set,  we  examine 
spreading  rates  between  20  and  120  km/my  and  we  vary  Q  by  changing  ATp  between 

100°C  and  200°C  (Table  2).  We  track  the  distribution  of  plume  material  by  introducing  a 
tracer  P  in  the  plume  and  using  our  tensor  diffusion  scheme  to  advect  P  passively  with  the 
mantle.  P=1  is  introduced  in  the  plume  source  column  to  represent  100%  plume  material, 
while  P=0  represents  0%  plume  material  and  100%  ambient  mantle.  We  define  Was  the 


along-axis  extent  to  which  the  depth-integrated  tracer  concentration 


1 


0.6D 


0.6  D 


jP(0,y,z)dz 


VO  j 

is  >0.05  (Fig.  2).  The  volume  flux  of  the  plume  is  measured  at  z=0.6D  by  integrating  the 
vertical  flow  of  the  plume  source  over  its  cross-sectional  area. 

For  calculations  that  include  thermal  buoyancy  only  without  latent  heat  loss  (set  A),  we 
find,  similar  to  ref.  [9,  1 1],  that  W depends  primarily  on  the  scaling  quantity  (Q/U)^^^,  and 


depends  secondarily  on  the  plume  buoyancy  number. 


Qp^aATp^ 
,  4Sri^U^  ^ 


as  defined  in  ref. 


[1 1],  and  on  the  ambient/plume  viscosity  ratio  r^rfo/rip,  at  z=0.5D.  A  modified  buoyancy 
number  which  depends  on  plume  viscosity  is  thus  (By).  The  best  fit  linear  regression 
function  obtained  by  fitting  linear  and  constant  coefficients  to  ln(By)  is 
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=  (11) 

Calculated  values  of  ^(2/(7)- 1/2  range  from  2.2  to  2.9  (Table  2)  with  a  mean  value  of 
2.50.  To  compare  our  results  directly  with  those  of  Ribe  et  al.  [11],  we  omit  the 
dependence  on  yand  incorporate  their  definition  of  Q  which  is  the  integrated  vertical  plume 
flux  weighted  by  plume  temperature  anomaly.  With  these  modifications  we  obtain  a  best- 
fit  linear  regression  of  W=2.%Q){QIU)^I'^B0  05  which  is  in  good  agreement  with  that  of  Ribe 
et  al.  [11]  of  While  the  scaling  and  exponential  factors  vary 

slightly  between  our  results  and  those  of  ref.  [9]  and  [11],  the  general  form  of  Eq.  (1 1)  is 
robust  and  insensitive  to  differences  in  far-field  experimental  boundary  conditions. 

For  calculations  of  thermal  buoyancy  with  latent  heat  loss  (set  B),  we  obtain  a  best-fit 
linear  regression  function, 

W  =  2.2l|^^J  (57)002.  (12) 

The  smaller  constant  and  exponential  coefficients  relative  to  those  in  Eq.  (11)  reflect  the 
inhibiting  effects  of  latent  heat  on  along-axis  plume  spreading.  The  average  values  of 
W(Q/U)-^/'2  for  experimental  set  B  is  2.29,  or  -92%  of  the  average  in  set  A. 

Addition  of  depletion  buoyancy  in  experimental  set  C  results  in  a  best-fit  regression 
function, 

W  =  2.  (57)0-04  (13) 

This  function  is  essentially  the  same  as  that  of  Eq.  (1 1)  for  set  A.  The  average  value  of 
W(Q/U)-^^'^  of  2.51  is  also  essentially  the  same  as  that  in  set  A.  The  further  addition  of 
melt  retention  (set  D)  does  not  change  this  relationship  significantly  as  shown  by  the 
similarity  in  regression  lines  of  set  C  and  set  D  (Fig.  4).  Thus,  the  effects  of  retention 
buoyancy  occurs  at  wavelengths  too  short  to  affect  the  full  width  W.  In  summary,  the 
effects  of  latent  heat  loss  to  inhibit  lateral  plume  spreading  are  approximately  balanced  by 
the  added  buoyancy  of  melt  depletion  which  enhances  plume  spreading. 

5.  Models  of  Iceland  and  the  Mid-Atlantic  Ridge 

We  next  investigate  models  of  mantle  flow  and  melting  beneath  Iceland,  a  relatively  well 
studied  example  of  a  ridge-centered  plume.  Our  objective  is  to  constrain  the  temperature 
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anomaly,  dimension,  and  volume  flux  of  the  Iceland  plume  by  comparing  theoretical  model 
predictions  with  observed  along-axis  variations  in  seismic  crustal  thickness,  topography, 
gravity,  and  basalt  geochemistry.  Previous  geophysical  studies  of  the  Iceland-MAR 
system  demonstrated  that  the  topographic  high  at  Iceland  coincides  with  a  low  in  mantle- 
Bouguer  gravity  anomaly  (MBA),  and  that  both  MBA  and  topographic  anomalies  can  be 
explained  by  the  combined  effects  of  anomalously  thick  crust  and  low  density  mantle 
generated  by  the  Iceland  plume  [27,  28],  MBA  are  calculated  by  subtracting  from  free-air 
gravity  the  attraction  of  seafloor  topography  and  the  crust-mantle  interface  assuming  a 
uniform  crustal  thickness  of  7  km  (e.g.,  [29,  30]).  Because  as  much  as  75%  of  the  along- 
axis  topographic  and  MBA  variations  may  arise  from  thickened  igneous  crust  [28,  31], 
crustal  thickness  calculations  are  an  important  link  between  our  models  and  surface 
observations. 

To  predict  crustal  thickness  from  mantle  melting  calculations,  we  assume  that  all  melt 
generated  within  200  km  of  the  ridge  axis  accretes  perpendicularly  to  the  ridge  axis  and  take 
the  top  of  our  numerical  box  to  be  the  isostatic  depth  of  the  seafloor  for  crust  of  normal 
thickness  (7  km).  The  crustal  thickness  as  a  function  of  along-axis  coordinate  y  is 
therefore 


Cr{y)  =  - 
U 


^^jjM(y)cb:dz. 


(14) 


We  take  the  top  of  our  model  to  be  the  isostatic  depth  of  the  seafloor  for  a  7-km-thick 
model  crust,  and  assume  isostatic  compensation  of  crustal  thickness  variations  that  deviate 
from  this  model  crust.  Consequently,  variations  in  crustal  thickness  impart  no  lithostatic 
pressure  variations  in  the  mantle.  To  prevent  melting  at  depths  shallower  than  the  isostatic 
base  of  the  thickened  Icelandic  crust  we  prohibit  melting  everywhere  at  depths  <28  km. 
Melting  may  stop  deeper,  however,  if  hydrothermal  cooling  is  important  [32]. 

To  calculate  isostatic  topography  of  the  seafloor,  we  consider  contributions  from  both 
the  crust  (Ahc)  and  mantle  (A/im).  In  calculating  Ahc,  we  assume  Airy  compensation  of  the 
crust  with  a  surface  density  contrast  of  (Pc-Pw)  for  the  submarine  portion  of  topography 
and  Pc  for  the  subaerial  portion.  The  crust  along  the  Reykjanes  and  Kolbeinsey  Ridges  is 
assumed  to  have  a  density  of  2800  kg/m^  except  within  500  km  of  the  plume  center,  where 
we  increase  it  linearly  to  a  maximum  of  3030  kg/m^  at  Iceland,  to  account  for  the  higher 
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MgO  content  of  the  Icelandic  crust  [33],  The  mantle  contribution  to  topography,  or 
dynamic  topography  is  calculated  from  vertical  normal  stress  at  the  top  layer  of  our  model. 


Ah^  = 


(Po  Pw^8 


(15) 


With  this  definition,  our  calculations  predict  seafloor  depths  to  increase  approximately  with 
the  square-root  of  distance  from  the  ridge-axis  which  is  consistent  with  lithospheric  half 
space  cooling  models  (e.g.  [34]).  In  addition  to  using  Ah^  to  predict  topography,  we  also 
use  A/zm  to  estimate  crustal  thickness  in  a  manner  independent  of  our  mantle  melting 
calculations.  This  "isostatic  crustal  thickness"  is  defined  as  the  isostatic  thickness  of  crust 
required  to  account  for  the  difference  between  the  observed  topography  and  Ahm- 

In  computing  MBA  we  again  consider  both  crustal  and  mantle  contributions.  The 
crustal  contribution  is  the  gravitational  signal  due  to  undulations  at  the  crust-mantle 


interface  that  deviate  from  the  constant  crustal  thickness  reference  model  originally  assumed 
in  generating  MBA.  For  these  calculations  we  employ  the  method  of  ref.  [35].  The  mantle 
contribution  to  gravity  is  calculated  by  integrating  the  contributions  from  lateral  density 
variations  at  each  model  layer  [29]. 


5.1  Broad  plume  source  model 

Our  first  model  of  the  Iceland-MAR  system,  much  like  that  of  Ribe  et  al.  [11], 
considers  a  broad  plume  source  with  a  relatively  small  temperature  anomaly  (model  Ice  1, 
a=300  km,  and  Arp=75°C)  rising  beneath  a  model  MAR  with  a  full-spreading  rate  of  19 
km/my  [36]  (Fig.  5a).  At  this  spreading  rate,  ro=1350°C  is  required  to  produce  a  ~7  km- 
thick,  normal  oceanic  crust.  The  calculation  that  includes  all  melting  effects  (model  Ice  Id) 
predicts  a  plume  volume  flux  of  -1.2x107  km^/my,  generating  an  along-axis  plume-head 
width  W of  -2300  km  (Fig.5a).  The  predicted  maximum  upwelling  rate  in  model  Ice  Id  is 
105  km/my,  which  is  >10  times  that  beneath  the  unaffected  portion  of  the  ridge  far  from  the 
plume.  The  predicted  upwelling  rate  averaged  through  the  melting  zone  in  plume  center  is 
20  km/my.  Melt  retention  buoyancy  contributes  minimal  effects  to  this  average  upwelling 
rate  and  thus  very  little  to  melting  rate. 

The  enhanced  upwelling  rate  in  the  plume  center,  combined  with  an  increase  in  total 
melt  extent  (23%  compared  to  13%  beneath  the  ridge  far  away  from  the  plume),  generates  a 
maximum  crustal  thickness  of  —30  km,  consistent  with  the  seismic  measurements  on 
Iceland  [37]  (Fig.  5b).  Along  the  length  of  the  Reykjanes  Ridge,  the  crustal  thickness 
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profiles  predicted  by  melting  in  model  Ice  Id  shows  striking  similarity  to  the  seismic 
measurements.  From  the  plume  center,  the  predicted  crust  first  thins  to  9.5  km  at  an  along- 
axis  distance  of  -300  km,  then  thickens  to  1 1  km  at  a  distance  of  -500  km,  and  finally 
tapers  to  a  thickness  of  6.7  km  at  a  distance  of  —1300  km.  The  predicted  local  minimum  in 
crustal  thickness  at  y-300  km  is  caused  by  a  reduced  mantle  upwelling  rate  at  the  plume 
edge  caused  by  the  rapid  vertical  flux  in  the  plume  center.  Melt  retention  does  not 
significantly  affect  crustal  thickness  because  the  predicted  0.5%  contrast  in  porosity 
between  the  plume  center  and  normal  sub-ridge  mantle  is  too  small  to  appreciably  enhance 
plume  upwelling  rate  in  the  shallowest  100  km,  where  melting  occurs.  The  isostatic  crustal 
thickness  profile  of  model  Ice  Id  also  shows  good  agreement  with  the  observed  crustal 
thickness  profile  (Fig.  5b).  The  excess  magmatic  flux  rate  required  to  sustain  the 
anomalous  (in  excess  of  a  7-km-thick  crust)  isostatic  crust  along  the  MAR,  1000  km  north 
and  south  of  Iceland,  is  2.33x10^  km^/my.  This  value  is  within  a  few  percent  of  the 
2.45x1 0^  km^/my  excess  crustal  production  rate  predicted  from  our  melting  model. 

The  predicted  topography  from  the  melting-model  crustal  profile  generates  70%  (-2.5 
km)  of  the  total  along-axis  topographic  anomaly  of  -3.5  km  (Fig.  5c).  We  predict  the 
remaining  30%  (-1  km)  of  topography  to  be  supported  by  dynamic  mantle  uplift  which  is 
obtained  with  a  P  value  of  0.024  [12,  17].  Of  mantle  dynamic  topography  Ah^,  thermal 
buoyancy  generates  -70%  while  depletion  and  retention  buoyancy  generate  the  remaining 
22  and  8%,  respectively.  The  predicted  total  amplitude  of  Ah^n  is  consistent  with  the  0.5- 
1.5  km  of  Eocene  uplift  as  inferred  from  sediment  core  analyses  [38]. 

The  mantle-Bouguer  anomaly  along  the  submarine  portions  of  the  ridge  is  also  matched 
well  by  predictions  of  model  Ice  Id  using  both  the  melting-model  and  isostatic  crust  (Fig. 
5d).  Similar  to  bathymetry,  the  crustal  MBA  accounts  for  most  (70%)  of  the  total  predicted 
anomaly  of  -330  mGal  with  the  mantle  contributing  the  remaining  30%.  Of  the  predicted 
mantle  gravity  signal,  75%  is  from  thermal  expansion,  while  20  and  5%  are  generated  by 
melt  depletion  and  retention,  respectively.  The  successful  predictions  of  both  topography 
and  MBA  support  the  hypothesis  that  these  anomalies  are  from  the  same  sources:  primarily 
crustal  thickness  variations  and  secondarily  density  variations  in  the  shallow  mantle. 

5.2  Narrow  plume  source  model 

Our  second  set  of  models  (Ice  2)  represent  another  end-member  possibility — that  of  a 
narrower  and  hotter  plume  source  (Fig.  6a;  <3=60  km,  and  Arp=170°C).  With  all  melting 
effects  included,  model  Ice  2d  predicts  a  plume  volume  flux  of  2.1x10^  km^/my  which 
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spreads  plume  material  to  a  full  width  W  of  870  km  along  the  ridge  axis.  The  maximum 
upwelling  rate  of  model  Ice  2d  is  283  km/my,  which  is  >2.5  times  greater  than  the 
maximum  upwelling  rate  in  the  broad  plume  source  (model  Ice  Id),  and  ~30  time  faster 
than  normal  ridge  upwelling  rates.  In  addition,  the  maximum  extent  of  melting  is  increased 
to  30%.  Thus  a  larger  volume  of  mantle  material  is  predicted  to  circulate  more  rapidly 
through  a  thicker  melting  zone  relative  to  that  of  Ice  Id,  which  results  in  melting  rates  an 
order  of  magnitude  greater  than  those  in  model  Ice  Id  (Fig.  6a).  For  the  model  without 
melt  retention  (model  Ice  2c),  the  melting-zone  averaged  upwelling  rate  is  63  km/my  and 
the  maximum  melting-model  crustal  thickness  is  147  km.  With  melt  retention  (model  Ice 
2d),  the  2.9%  porosity  in  the  plume  is  sufficient  to  increase  the  predicted  melting-zone- 
averaged  upwelling  rate  to  80  km/my  and  the  maximum  melting-model  crustal  thickness  to 
166  km  (Fig  6b).  In  model  Ice  2d,  the  melting-model  crust  thins  to  3  km  at  an  along-axis 
distance  of  120  km,  where  upwelling  and  thus  melting  rate  is  strongly  reduced  at  the  edge 
of  the  rapidly  upwelling  plume  center  (Fig.  6a). 

The  high  maximum  crustal  thicknesses  predicted  by  the  narrow  plume  source,  melting 
model  drastically  exceed  calculations  of  previous  studies  that  assumed  passive  mantle 
upwelling  (e.g.  ref.  [28,  39])  and  drastically  exceed  the  observed  crustal  thicknesses  (Fig 
6b).  The  resulting  topographic  and  MBA  anomalies  also  fail  to  match  the  observations 
(Fig.  6c,  d).  The  isostatic  crustal  profile,  on  the  other  hand,  yields  predictions  in  much 
better  agreement  with  the  observed  crustal  thicknesses  (Fig  6b),  topography  (by  definition) 
(Fig.  6c),  and  MBA  (Fig.  6d)  along  the  ridge  axis.  Thus,  if  the  Iceland  plume  is 
comparable  in  radius  and  temperature  to  our  narrow  plume  source  model,  a  substantial 
portion  of  the  melt  produced  beneath  Iceland  must  accrete  more  uniformly  along-axis  than 
our  melting-model  crust,  much  like  our  isostatic  crustal  profile.  This  condition  suggests 
melt  migration  and/or  lower  crustal  ductile  flow  [40]  occurs  over  distances  of  several 
hundreds  of  kilometers  away  from  Iceland  along  the  Reykjanes  and  Kolbeinsey  Ridges. 

Because  the  mechanisms  of  along  ridge-axis  melt  transport  are  poorly  understood,  we 
do  not  attempt  to  model  this  process  in  this  study.  Instead,  we  assume  a  priori  that  along- 
axis  melt  redistribution  does  occur  and  that  the  end  result  of  this  process  leads  to  the 
isostatic  crustal  profile.  In  arriving  at  our  final  Ice  2  models,  we  thus  sought  values  of  ATp 
and  a  such  that  the  total  volume  rate  of  melt  produced  by  the  melting  model  matched  that 
required  to  sustain  the  isostatic  crustal  profile.  The  best  solutions  of  AT'p=170°C  and  a=60 
km  yield  a  total  excess  melt  production  rate  of  2.54x10^  km^/my  (model  Ice  2d),  which  is 
within  1  %  of  that  required  of  the  isostatic  crustal  profile. 
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In  these  narrower,  hotter  plume  source  models,  the  mantle  contribution  to  topography 
and  gravity  relative  to  the  crustal  contribution  becomes  much  larger  than  in  the  broader, 
cooler  source  models.  For  example,  model  Ice  2d  predicts  a  mantle  topographic  uplift  that 
is  51%  (1.8  km)  of  the  observed  along-axis  topographic  anomaly  (Fig.  6c),  and  a  mantle 
contribution  to  MBA  that  is  48%  (158  mGal)  of  the  observed  MBA  variation  (Fig.  6d). 
The  crust  therefore  generates  only  49  and  52%  of  the  total  topographic  and  MBA 
variations,  respectively.  Calculations  also  predict  the  importance  of  melt-related  buoyancy 
to  the  mantle  anomalies  to  be  significantly  greater  for  these  hotter  plume  source  models 
relative  to  the  cooler  source  models.  Thermal  buoyancy  is  predicted  to  produce  47%  of 
A/jm,  and  60%  of  the  mantle  MBA  variation;  melt  depletion  produces  39%  of  A/zm  and 
25%  of  the  mantle  MBA;  and  melt  retention  produces  the  remaining  14%  of  Ah^  and  15% 
of  the  mantle  MBA  variation. 

5.3  Reykjanes  Ridge  bathymetric  swell 

Similar  to  along-axis  topography,  we  predict  map  view  topography  by  adding  mantle 
dynamic  topography  (Eq.  (15))  and  isostatic  topography  of  the  crust  considering  only 
along-axis  variations  in  crustal  thickness.  For  model  Ice  Id,  we  use  the  melting-model 
crust  and  for  model  Ice  2d,  we  use  the  isostatic  crust.  Fig.  7  illustrates  the  observed 
topography  in  map  view  along  the  Reykjanes  Ridge  south  of  Iceland,  and  predictions  of 
models  Ice  Id  and  Ice  2d.  The  similarity  between  the  predictions  and  observations  at  broad 
wavelengths  (>~500  km)  are  compelling:  both  models  predict  the  ~2.0  km  across-axis 
decrease  in  broad  wavelength  topography  between  Iceland  and  an  across-axis  distance  of 
400  km  away  from  the  ridge-axis,  and  both  predict  the  south-pointing  swell,  elongated 
along  the  Reykjanes  Ridge.  As  demonstrated  above,  the  southward  deepening  of  the  ridge 
axis  reflects  crustal  thinning  and  mantle  density  increase  with  distance  from  the  Iceland 
plume  source.  But  perpendicular  to  the  ridge-axis,  seafloor  topography  is  dominated  by 
the  subsidence  of  the  cooling  lithosphere.  Thus,  contrary  to  previous  notions  (e.g.  [5,  6]), 
the  regional  bathymetric  swell  does  not  require  a  pipe-like  flow  of  plume  material  along  the 
ridge  axis.  Instead,  we  predict  the  plume  head  to  spread  radially  and  explain  the  general 
shape  of  the  elongated  Icelandic  swell  as  the  superposition  of  radial  plume  spreading  and 
across-axis  lithospheric  cooling.  The  models  presented  in  this  study,  however,  do  not 
consider  time-dependent  variations  in  crustal  accretion  which  may  also  contribute  to  across- 
axis  topographic  variations. 


56 


5.4  Rare-earth  element  and  isotopic  anomalies 

A  potentially  useful  independent  constraint  on  melting  depth  and  extents,  which  reflect 
mantle  temperature,  is  rare-earth  element  (REE)  concentrations  of  axial  basalts.  A  simple 
comparison  can  be  made  with  previous  inversions  of  melt  fraction  versus  depth  as 
calculated  by  White  et  al.  [33].  At  the  plume  center,  our  broad  plume  source  model  (Ice 
Id)  and  narrow  plume  source  model  (Ice  2d)  predict  melt  fractions  that  are  lower  and 
higher,  respectively,  than  White  et  al.'s  [33]  inversions  for  Krafla  volcano  on  Iceland  (Fig. 
8a).  The  potential  temperature  of  the  Iceland  plume-source,  therefore,  is  likely  to  be  1425- 
1520°C  as  represented  by  our  two  end-member  models.  At  -550  km  from  Iceland  on  the 
Reykjanes  Ridge,  model  Ice  Id  predicts  melting  depths  and  extents  closely  matching  those 
obtained  from  the  REE  inversions  [33]  (Fig.  8b).  Model  Ice  2d,  however,  underpredicts 
the  extents  and  depths  because  plume  material  from  our  narrow  plume  source  did  not 
spread  to  this  along-axis  distance.  Thus,  in  order  to  explain  the  REE  composition  of 
basalts  sampled  550  km  away  from  Iceland,  once  again  our  model  Ice  2d  seems  to  require 
plume-derived  melts  to  migrate  substantially  along  the  Reykjanes  Ridge  axis. 

While  REE  concentrations  reflect  melting  process  beneath  Iceland,  Sr  isotope  ratios 
may  reflect  the  concentration  of  the  plume  source  material  relative  to  that  of  normal  mid¬ 
ocean  ridge  basalts  (MORE).  Schilling  [8, 41]  interprets  the  peak  in  ^^Sr/^^Sr  at  Iceland  to 
mark  the  center  of  the  Iceland  plume,  where  the  plume  source  concentration  is  highest,  and 
interprets  the  decrease  in  ^^Sr/^^Sr  north  and  south  of  Iceland  to  reflect  a  decrease  in 
percent  of  plume  material  comprising  the  mantle  melt  source. 

To  address  questions  of  where  and  how  plume— MORE  mixing  occurs,  we  calculate  the 
fraction  of  plume  tracer  P  in  accumulated  melts  along  the  model  ridge  axis  (neglecting 
along-axis  melt  migration)  (Fig  8c).  At  each  numerical  grid  where  new  melt  is  generated,  P 
is  weighted  by  melting  rate.  We  then  integrated  over  each  ridge-perpendicular  plane  to 
compute  a  weighted  mean  value  ( F )  for  each  point  along  the  ridge  axis, 

_  \P{x,y,z)M{x,y,z)AxAz 

P{y)  =  ^ - r-. - ■  (16) 

I  M{x,y,z)dx&z 

This  calculation  thus  approximates  the  plume  concentration  of  pooled  melts  along  the  ridge 
axis.  For  example,  F =1.0  indicates  that  all  of  the  melt  generated  in  a  plane  perpendicular 
to  that  point  of  the  ridge  is  entirely  plume-source  derived.  Likewise,  F  =0.0  indicates  that 
none  of  the  melts  are  plume  derived  and  0.0<  F  <1.0  indicates  plume-MORE  mixing. 
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Model  Ice  Id  predicts  an  along-axis  geochemical  plume  width  of  >2000  km, 
significantly  greater  than  that  suggested  by  the  87Sr/86Sr  anomaly.  Ice  2d  on  the  other 
hand  predicts  a  width  of  -1000  km  which  is  more  consistent  with  that  of  the  ^^Sr/^^Sr 
anomaly;  however,  its  profile  in  P  would  likely  be  broader  if  along-axis  melt  migration 
were  considered.  Both  model  Ice  Id  and  Ice  2d  predict  that  the  melts  are  entirely  plume 
derived  ( P =1 .0)  over  most  of  the  plume  width,  and  become  fully  ambient  mantle  derived 
( P  =0.0)  within  200-300  km  of  the  edge  of  the  plume.  These  results  suggest  that  within 
most  of  the  plume  affected  portion  of  the  ridge,  very  little  mixing  occurs  between  plume 
and  ambient  source  material  in  the  shallow  mantle.  Thus,  if  the  gradients  in  ^^Sr/^^Sr 
away  from  Iceland  reflect  plume-MORB  mixing,  it  most  likely  occurs  deeper  in  the  mantle, 
possibly  by  ambient  mantle  entrainment  of  the  ascending  plume  (e.g.  [42]). 

5.5  Predictions  of  P-wave  seismic  velocity  anomalies 

Observations  of  compressional  wave  (P-wave)  seismic  travel  time  variations  and 
associated  mantle  P-wave  velocity  variations  provide  critical  constraints  on  mantle 
properties  beneath  Iceland.  To  predict  P-wave  seismic  velocity  anomalies,  we  assume  a 
reference  P-wave  velocity  of  8  km/s,  which  decreases  by  6.25x10'^%  for  each  1°C  increase 
in  mantle  temperature,  increases  by  0.1%  for  each  1%  increase  in  depletion,  and  decreases 
by  1.25%  for  each  1%  increase  in  pore  volume  [43].  We  also  predict  P-wave  travel-time 
residuals  by  calculating  travel  times  of  seismic  rays  passing  vertically  through  the  400  km 
thickness  of  our  mantle  models. 

The  broad  plume  source  model  (Ice  Id)  predicts  a  maximum  decrease  in  P-wave 
velocity  below  the  melting  region  of  -0.5%  relative  to  the  surrounding  mantle.  In  the 
melting  region,  the  predicted  P-wave  velocity  anomaly  diminishes  because  the  velocity¬ 
enhancing  effects  of  latent  heat  loss  and  melt  depletion  exceed  the  velocity-reducing  effect 
of  melt  retention  (Fig.  9a).  The  corresponding  travel-time  delay  for  vertically  passing  rays 
is  predicted  to  be  +0.23  s  at  the  plume  center  and  decrease  to  zero  at  an  along-axis  distance 
of  -1200  km.  The  contributions  to  travel-time  delay  above  the  plume  center  are  +0.25  s 
from  excess  mantle  temperature,  -0.09  s  from  melt  depletion,  and  +0.07  s  from  melt 
retention.  Across  the  ridge-axis,  lithospheric  cooling  dominates,  resulting  in  a  predicted 
travel-time  difference  of  0.5  s  between  the  plume  center  and  at  an  across-axis  distance  of 
400  km.  The  broad  plume  source  model  thus  predicts  only  a  gradual  decrease  in  travel¬ 
time  delay  across  the  ridge  axis  and  even  smaller  variations  along  the  ridge  axis. 
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In  contrast,  the  narrow  plume  source  of  model  Ice  2d  predicts  significantly  larger 
amplitudes  of  P-wave  anomalies  over  a  much  narrower  lateral  extent.  Below  the  melting 
zone,  the  170°C  plume  temperature  anomaly  reduces  calculated  P-wave  velocities  by  more 
than  1%.  In  the  melt  zone,  however,  the  P-wave  velocities  are  reduced  to  as  much  as  2% 
due  to  the  2.9%  melt  retention  (Fig.  9b).  Along  the  ridge  axis,  the  travel-time  delay  for 
vertically  passing  rays  is  predicted  to  be  -1-0.75  s  at  the  plume  center  and  to  decrease  by 
0.85  s  within  -80  km.  Approximately  half  of  this  travel-time  residual  is  predicted  to  arise 
in  the  high-porosity  melt  zone  in  the  shallow  mantle.  Across  the  ridge  axis,  the  additional 
effect  of  lithospheric  cooling  yields  a  predicted  travel-time  difference  of  1  s  within  -80  km 
of  the  plume  center  and  a  travel-time  difference  of  1 .2  s  over  an  across-axis  distance  of  400 
km.  Preliminary  results  of  the  ongoing  ICEMELT  experiment  at  Iceland  have  revealed 
azimuthal  variations  in  P-wave  travel  times  as  high  as  1  s  within  100  km  of  the  ridge  axis 
[44],  suggesting  that  the  narrow  plume  source  model  better  represents  Iceland  than  does  the 
broad  plume  source  model. 

6.  Discussion 

6.1  Importance  of  melting  effects 

The  importance  of  melting  effects  on  mantle  flow,  melt  production,  and  surface 
observables  are  summarized  in  Fig.  10.  Mantle  melting  generates  appreciable  effects  on 
mantle  properties;  however,  over  the  range  of  plume  viscosities  considered  in  our  models, 
the  effect  of  latent  heat  loss  on  mantle  flow  largely  cancels  the  effects  of  depletion  and 
retention  buoyancy.  As  a  result,  the  combined  effects  of  these  factors  on  mantle  flow  are 
small  as  reflected  in  the  small  changes  in  the  predicted  values  of  W{Q/U)-^I'^  (Figs.  4  and 
10).  Similarly,  when  plume  temperature  anomalies  are  mild  as  in  the  Ice  1  models,  the 
melting-related  factors  have  only  second  order  effects  on  upwelling  rate  as  reflected  in 
small  changes  in  the  predicted  crustal  thickness  (Fig.  10).  When  plume  temperature 
anomalies  are  larger,  however,  as  in  the  Ice  2  models,  melt  retention  may  enhance  the 
predicted  crustal  thickness  by  20%  relative  to  calculations  that  do  not  include  retention. 

Contrasting  with  their  mild  influence  on  mantle  flow,  the  melting-related  factors  have 
substantial  effects  on  the  predicted  geophysical  observables  and  these  effects  increase  with 
increasing  plume  temperature  (Fig.  10).  For  mantle  contributions  to  topography  and  MBA, 
latent  heat  loss  reduces  the  amplitudes  of  predicted  anomalies  by  20-40%  relative  to 
calculations  without  latent  heat  loss.  Depletion  buoyancy  increases  predicted  mantle 
topographic  anomalies  and  MBA  by  10-65%  relative  to  calculations  without  depletion. 
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while  retention  buoyancy  increases  predicted  anomalies  by  5-25%  relative  to  calculations 
without  retention.  Melting  effects  on  P-wave  delay-time  are  also  important:  Latent  heat 
loss  decreases  predicted  delay-time  by  -13%,  melt  depletion  decreases  delay-time  by  20- 
30%,  but  melt  retention  increases  delay-time  by  20-60%.  It  is  thus  important  to  consider 
melting  effects  on  mantle  properties  when  predicting  geophysical  observables. 

6.2  Model  uncertanties 

Because  melting-related  factors  do  not  affect  significantly  large-scale  mantle  flow, 
uncertainties  associated  with  our  melt  calculations  such  as  the  assumed  batch  melting  [25], 
our  choice  of  P  values,  and  the  melt  porosity  calculations,  are  likely  to  have  only  secondary 
effect  on  our  estimates  of  plume  source  radius  and  temperature.  By  far  the  most  important 
uncertainty  in  this  regard  is  mantle  rheology.  The  reference  mantle  viscosity  rjo  controls 
directly  the  rate  of  mantle  upwelling  in  response  to  density  variations  (Eq.  (9)).  But 
unfortunately,  viscosity  beneath  ridges  is  not  known  to  within  one  or  even  two  orders  of 
magnitude  (e.g.  [45]).  One  mechanism  that  may  yield  a  substantially  higher  viscosity  than 
that  we  have  assumed  is  dehydration  at  the  onset  of  melting  [45].  A  higher  melting  zone 
viscosity,  for  example,  would  most  likely  require  a  greater  temperature  anomaly  of  the 
broad  plume  source  model  to  explain  the  geophysical  observations,  or  require  a  greater 
source  radius  and  less  along-axis  melt  redistribution  of  our  narrow  plume  source  model  to 
explain  the  observations.  Thus,  because  of  the  uncertainty  of  viscosity,  our  Iceland  plume 
models  are  not  unique.  However,  they  do  provide  reasonable  bounds  on  the  plume  source 
radius  and  temperature  given  the  similarities  between  model  predictions  and  the  variety  of 
geophysical  and  geochemical  observations  considered. 

6.3  Plume  volume  flux  estimates 

Still,  it  may  be  possible  to  constraint  plume  volume  flux  independent  of  ambient 
viscosity  based  on  the  observed  MBA  and  bathymetric  anomaly  widths  and  the  theoretical 
relationship  between  flux  and  W  (i.e.  Eq.  (13)).  The  use  of  Eq.  (13)  to  infer  plume 
volume  flux  is  valid  if  the  surface  anomaly  widths  reflects  directly  the  along-axis  plume 
width  in  the  mantle,  which  would  be  the  case  if  along-axis  melt  migration  is  negligible  as 
assumed  in  the  Ice  I  models.  The  flux  required  to  match  the  along-axis  MBA  and 
bathymetric  anomaly  widths  as  predicted  from  model  Ice  Id  is  1.2x107  km^/my.  This 
flux,  however,  is  several  times  larger  than  previous  estimates  of  the  Iceland  plume  of 
2x106  km/my  [46],  1.43x106  km/my  [8],  and  2.2x106  km/my  [28].  If  on  the  other  hand. 
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along-axis  melt  migration  is  important  as  suggested  for  the  Ice  2  models,  we  can  not  use 
Eq.  (13)  to  constrain  the  Iceland  plume  volume  flux  independent  of  rjo-  We  must  therefore 
rely  on  the  fact  that  our  melt  production  rate  estimates  are  consistent  with  the  total  volume 
of  observed  excess  crust  as  we  did  for  the  Ice  2  models.  Indeed,  model  Ice  2d  predicts  a 
plume  volume  flux  of  2.1x106  km^/my  which  is  more  consistent  with  the  above  estimates 
of  the  previous  studies.  An  intriguing  new  question  arising  from  this  narrow  plume  source 
model  is,  what  specific  mechanisms  may  allow  melt  generated  beneath  Iceland  to  migrate 
hundreds  of  kilometers  along-axis?  Possible  evidence  for  such  melt  transport  may  include 
the  V-shaped  axial  bathymetric  highs  propagating  away  from  Iceland  along  the  Reykjanes 
Ridge  as  first  noted  by  Vogt  [47]  in  1971. 

7.  Conclusions 

We  have  investigated  the  dynamics  of  mantle  flow  and  melting  of  a  ridge-centered 
plume  using  three-dimensional,  variable-viscosity  models  with  focus  on  three  buoyancy 
sources:  temperature,  melt  depletion,  and  melt  retention.  When  all  melting  effects  are 
considered,  the  relationship  between  along-axis  plume  width  W,  plume  volume  flux  Q,  full 
spreading  rate  U,  buoyancy  number  B,  and  ambient/plume  viscosity  ratio  7,  is  best 
parameterized  by  Calculations  that  include  melting  yield  a 

similar  relationship  to  those  that  do  not  include  melting  because  of  the  competing  effects  of 
latent  heat  loss  and  depletion  buoyancy.  We  propose  two  end-member  models  for  the 
Iceland  plume  beneath  the  MAR.  The  broad  plume  source  of  radius=300  km  represents  a 
low  temperature  (Arp=75°C)  and  high  flux  (2=1-2x10'^  km^/my)  end-member,  while  the 
narrow  plume  source  of  radius=60  km  represents  a  high  temperature  (Arp=170°C)  and 
low  flux  (2=2.1x106  km^/my)  end-member.  The  broad  plume  source  predicts 
successfully  the  observed  along-axis  variations  in  seismic  crustal  thickness,  topography, 
and  mantle-Bouguer  gravity  anomalies;  whereas  the  narrow  source  model  predicts 
adequately  the  total  excess  crustal  production  rate  (2.5x106  km^/my)  but  requires  extensive 
melt  migration  and/or  lower  crustal  ductile  flow  to  occur  over  hundreds  of  kilometers  along 
the  MAR  in  order  to  explain  the  geophysical  and  geochemical  observations.  Our 
calculations  predict  that  plume  spreading  away  from  the  plume  center  is  radially  symmetric 
rather  than  channeled  preferentially  along  the  ridge  axis.  The  elongated  bathymetric  swell 
along  the  Reykjanes  Ridge  can  be  explained  by  rapid  off-axis  subsidence  due  to 
lithospheric  cooling  superimposed  on  a  broader  hotspot  swell.  Both  the  broad  and  narrow 
plume  source  models  predict  very  little  mixing  between  the  plume  and  MORB  sources  in 


61 


the  shallow  mantle;  hence,  we  suggest  that  mixing  may  occur  deeper  in  the  mantle  possibly 
due  to  entrainment  of  the  isotopically  depleted  portion  of  the  mantle  by  the  rising  mantle 
plume.  Our  narrow  plume  source  model  predicts  seismic  P-wave  velocity  variations  more 
consistent  with  recent  seismic  observations  beneath  Iceland,  suggesting  that  this  model  may 
better  represent  the  Iceland  plume. 
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Table  1.  Notation 


Variable 

Meanine 

Value 

Units 

a 

plume  radius 

km 

b 

grain  size 

3x10“ 

m 

B 

buoyancy  number 

specific  heat 

1000 

J  kg-'  °C-' 

D 

fluid  depth 

400 

km 

E 

activation  energy 

1.9x105 

J 

8 

acceleration  of  gravity 

9.8 

m/s 

Ah^ 

isostatic  crustal  topography 

km 

Ah^ 

mantle  dynamic  topography 

km 

K 

mantle  permeability 

m^ 

M 

melt  fraction 

wt% 

P 

pressure 

Pa 

P 

plume  tracer  concentration 

Q 

volumetric  plume  flux 

kmVmy 

R 

gas  constant 

8.314 

J  K-'  mol-' 

Ra 

thermal  Rayleigh  number 

Ra^ 

depletion  Rayleigh  number 

Ra^ 

retention  Rayleigh  number 

AS 

entropy  change  on  melting 

400 

J  kg-'  °C 

T 

mantle  potential  temperature 

°C 

Tr 

mantle  real  temperature 

K 

AT^ 

plume  temperature  anomaly 

U 

mantle  flow  rate  vector 

km/my' 

U 

ridge  full  spreading  rate 

km/my 

V 

activation  volume 

4x10-6 

m^ 

w 

along-axis  plume  width 

km 

X 

melt  depletion 

wt% 

a 

coefficient  of  thermal  expansion 

3.4x10-5 

K-» 

p 

coefficient  of  depletion  density  reduction 

r 

rio'rip 

K 

thermal  diffusivity 

31 

km^/my 

rt 

viscosity 

Pa  s 

Ro 

reference  viscosity 

Pa  s 

rip 

plume  viscosity  at  0.5D 

Pa  s 

Vm 

melt  viscosity 

1.0 

Pa  s 

CO 

vertical  melt  flow  rate 

km/my 

p 

mantle  density 

kg/m^ 

Pc 

crust  density 

2800^3030 

kg/m^ 

Pm 

melt  density 

2900 

kg/m^ 

Po 

mantle  reference  density 

3300 

kg/m^ 

Pw 

water  density 

1000 

kg/m^ 
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Table  2a.  Parameters  and  results  of  experimental  set  A:  thermal  buoyancy  without  latent  heat  loss 


Model 

grid  size,  dx-dy  (km) 

U  (km/my) 

B 

r 

Q  (lO^km-^/my) 

W 

w{Q/ur'^ 

1  a 

12.5-12.5 

20 

100 

1.729 

2.352 

0.974 

512 

2.322 

2a 

12.5-6.25 

120 

100 

0.0522 

2.352 

1.059 

219 

2.328 

3a 

12.5-6.25 

60 

100 

0.195 

2.352 

0.987 

281 

2.193 

4a 

12.5-6.25 

40 

100 

0.460 

2.352 

1.038 

362 

2.251 

5a 

12.5-12.5 

20 

200 

6.977 

5.054 

1.965 

938 

2.991 

6a 

12.5-6.25 

120 

200 

0.244 

5.054 

2.478 

331 

2.305 

7a 

12.5-6.25 

60 

200 

0.746 

5.054 

1.892 

419 

2.358 

8a 

12.5-12.5 

50 

200 

1.468 

5.054 

2.585 

662 

2.914 

Ice  la 

12.5-12.5 

19 

75 

36.553 

1.849 

12.39 

2312 

2.864 

Ice  2a 

10.0-10.0 

19 

170 

14.866 

3.757 

2.186 

830 

2.447 

For  all  experiments  vertical  grid  separation  dz  is  8  km. 

Models  1-8:  /?=0.06,  Ar=1300°C,  s,  Ra=0,9\5  x  \0^,Rax=\35  x  10^ /?a<^=2.75  x  10^  a=70 

km. 


Models  Ice  1-2:  ^=0.024,  Ar=I350°C,  77o^5xlO*^Pa  s,  Ra=\.90  x  10^  Rax=\A2  x  10^  Ra^SJO  x  10^ 
Model  Ice  1:  a=300  km,  Arp=:75°C. 

Model  Ice  2:  a=60  km,  Arp=170°C. 


Table  2b.  Parameters  and  results  of  experimental  set  B:  thermal  buoyancy  with  latent  heat  loss 


Model 

grid  size,  dx-dy  (km) 

U  (km/my) 

B 

r 

Q  (10^  km^/my) 

W^Q/U)-''^ 

lb 

12.5-12.5 

20 

100 

1,720 

2.352 

0.969 

488 

2.215 

2b 

12.5-6.25 

120 

100 

0.0463 

2.352 

0.939 

206 

2.331 

3b 

12.5-6.25 

60 

100 

0.166 

2.352 

0.843 

256 

2.162 

4b 

12.5-12.5 

40 

100 

0.389 

2.352 

0.876 

312 

2.111 

5b 

12.5-12.5 

20 

200 

7.001 

5.054 

1.972 

838 

2.667 

6b 

12.5-6.25 

120 

200 

0.203 

5.054 

2.059 

281 

2.147 

7b 

12.5-6.25 

60 

200 

0.750 

5.054 

1.901 

369 

2.072 

8b 

12.5-12.5 

50 

200 

1.193 

5.054 

2.100 

462 

2.257 

Ice  lb 

12.5-12.5 

19 

75 

34.774 

1.849 

1 1.79 

2212 

2.809 

Ice  2b 

10. 0-10.0 

19 

170 

14.857 

3.757 

2.185 

710 

2.094 
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Table  2c.  Parameters  and  results  of  experimental  set  C:  thermal  +  depletion  buoyancy  with  latent  heat  loss 


Model 

grid  size,  dx-dy  (km) 

U  (km/my) 

B 

/ 

Q  (lO^km-^my) 

W{Q/U)-‘^ 

Ic 

12.5-12.5 

20 

100 

1.427 

2.352 

0.804 

488 

2.431 

2c 

12.5-6.25 

120 

100 

0.0452 

2.352 

0.917 

206 

2.360 

3c 

12.5-6.25 

60 

100 

0.164 

2.352 

0.834 

256 

2.173 

4c 

12.5-12.5 

40 

100 

0.385 

2.352 

0.868 

338 

2.291 

5c 

12.5-12.5 

20 

200 

6.788 

5.054 

1.912 

988 

3.194 

6c 

12.5-6.25 

120 

200 

0.188 

5.054 

1.902 

281 

2.234 

7c 

12.5-6.25 

60 

200 

0.724 

5.054 

1.835 

406 

2.323 

8c 

12.5-12.5 

50 

200 

1.159 

5.054 

2.040 

538 

2.661 

Ice  Ic 

12.5-12.5 

19 

75 

34.1 17 

1.849 

1 1.56 

2288 

2.932 

Ice  2c 

lO.O-IO.O 

19 

170 

14.597 

3.757 

2.147 

830 

2.469 

Table  2d. 

Parameters  and  results  of  experimental  set  D 

:  thermal  +  depletion  +  retention  buoyancy  with  latent  heat 
loss 

Model 

grid  size,  dx-dy  (km) 

U  (km/my) 

B 

7 

g  (10^  km'^my) 

W 

WiQ/U)-‘^ 

id 

12.5-12.5 

20 

100 

1.433 

2.352 

0.808 

488 

2.426 

2d 

12.5-6.25 

120 

100 

0.0452 

2.352 

0.917 

206 

2.359 

3d 

12.5-6.25 

60 

100 

0.1648 

2.352 

0.836 

256 

2.171 

4d 

12.5-12.5 

40 

100 

0.358 

2.352 

0.868 

338 

2.292 

5d 

12.5-12.5 

20 

200 

6.764 

5.054 

1.905 

913 

2.956 

6d 

12.5-6.25 

120 

200 

0.193 

5.054 

1.961 

281 

2.200 

7d 

12.5-6.25 

60 

200 

0.718 

5.054 

1.820 

419 

2.404 

8d 

12.5-12.5 

50 

200 

1.139 

5.054 

2.005 

563 

2.809 

Ice  Id 

12.5-12.5 

19 

75 

34.089 

1.849 

11.56 

2288 

2.933 

Ice  2d 

10.0-10.0 

19 

170 

14.501 

3.757 

2.133 

870 

2.597 
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Fig.  1.  Combined  shipboard  and  Etopo5  bathymetry  map  (contour  interval  of  0.5  km) 
showing  Iceland  (65°N,  18°W)  and  the  Reykjanes  (south  of  Iceland)  and  Kolbeinsey  (north 
of  Iceland)  Ridges.  Bold  lines  marks  the  ridge  axes.  This  figure  and  Figs.  4,  7,  8,  and  10 
were  produced  using  the  GMT  software  package  [48]. 
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Fig.  2.  Perspective  diagram  illustrating  steady-state  flow  (small  arrows)  and  potential 
temperature  (shaded  and  contoured  at  100°C-intervals)  fields  of  an  example  calculation  that 
considers  thermal  buoyancy  only  and  no  melting  effects  (model  5a).  Vertical  plane  on  the 
right  is  a  depth  cross-section  along  the  ridge  axis  (jc=0),  while  the  vertical  plane  to  the  left 
is  a  depth  cross-section  perpendicular  to  the  ridge  axis  (y=0).  Top  plot  shows  depth- 
averaged  plume  tracer  concentration  P  along  the  ridge  axis  which  we  used  to  define  plume 
width  W.  Both  top  (z=0)  and  bottom  (z=D)  boundaries  are  isothermal  planes  with  the 
bottom,  a  free  slip  boundary  and  the  top,  fixed  at  a  horizontal  velocity  of  0.5(7  (large 
horizontal  arrow).  All  boundaries  are  closed  to  flow  both  in  and  out  of  the  numerical  box, 
thus  material  flows  downward  at  the  end  of  the  box  opposite  the  ridge  (x=800  km)  and 
recirculates  toward  the  ridge  axis  along  the  base  of  the  box.  The  effect  of  this  recirculation 
on  the  interaction  between  plume  and  ridge  are  insignificant. 
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Fig.  3.  Perspective  views  of  depth  cross-sections  showing  %  density  reduction  in  the 
mantle  due  to  (a)  thermal  buoyancy  with  latent  heat  loss  {ATp=200°C)  (model  5b),  (b)  plus 
melt  depletion  buoyancy  (model  5c),  and  (c)  plus  melt  retention  buoyancy  (model  5d). 
Contour  interval  is  0.5%.  Vectors  in  (a)  show  mantle  flow.  Vectors  in  (b)  show  the 
differences  between  flows  with  and  without  melt  depletion  buoyancy.  Vectors  in  (c)  show 
the  difference  between  flows  with  and  without  melt  retention  buoyancy.  Downward 
pointing  vectors  in  (b)  and  (c)  illustrate  reduced  upwelling,  not  downwelling. 


Fig.  4.  Numerical  results  (dots)  of  calculations  with  all  melting  effects  included  (set  D). 
The  two  Iceland  models  are  circled.  The  solid  black  line  is  the  best-fit  linear  regression 
shown  by  Eq.  (13)  which  yields  a  standard  deviation  misfit  that  is  7%  of  the  median  value 
of  W{QIU)'^l'^.  Also  shown  are  corresponding  linear  regressions  of  calculations  of  thermal 
buoyancy  without  latent  heat  loss  (set  A,  gray),  thermal  buoyancy  with  latent  heat  loss  (set 
B,  dotted),  and  additional  buoyancy  from  melt  depletion  (set  C,  dashed). 


Figure  4 


77 


Fig.  5.  (a)  Perspective  diagram  of  model  Ice  Id  (broad  plume  source)  shaded  according  to 
temperature.  Black  contours  are  depletion  (contour  interval  is  5%)  and  white  contours  are 
melting  rates  of  0.01,  0.03,  and  0.05  my^.  (b)  Comparison  between  model  Ice  Id 
melting-model  crust  (solid)  and  isostatic  crust  (dashed),  and  seismic  crustal  thickness 
measurements  along  the  Reykjanes  Ridge  (dots)  and  at  older  seafloor  near  the  continental 
margins  (triangles)  from  ref.  [37].  (c  and  d)  Comparison  between  the  observed  bathymetry 
(thick  gray  curve  in  c)  and  MBA  (thick  gray  curve  in  d)  along  the  MAR  and  predicted 
profiles  of  model  Ice  Id  using  the  melting-model  crust  (bold  curves  in  c  and  d)  and  isostatic 
crust  (thick  dashed  curved  in  d).  Also  shown  are  predicted  mantle  components  due  to 
various  mantle  density  sources  as  labeled.  Bathymetry  data  and  MBA  are  from  ref  [28]. 
We  do  not  consider  on-land  gravity  of  Iceland. 
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Fig.  6.  Same  as  Fig.  5  but  for  Ice  2  models  (narrow  plume  source).  Symbols  are  the  same 
as  in  Fig.  5  except  melting  rate  contours  in  (a)  are  0.01,  0.03,  0.05,  0.2,  0.4  myF 
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Fig.  7.  (a)  Observed  topography  of  Iceland  and  the  Reykjanes  Ridge  (oblique  Mercator 
projection),  (b)  Mantle  +  crustal  topography  predicted  from  our  broad  plume  source  model 
Ice  1  d  using  the  melting-model  crust,  (c)  Mantle  -l-  crustal  topography  predicted  from  our 
narrow  plume  source  model  Ice  2d  using  the  isostatic  crust. 
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a)  Observed  Topography  b)  Model  Ice  1d 
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Fig.  8.  (a  and  b)  Comparison  between  White  et  al.’s  [33]  REE  inversion  of  melt  fraction 
(gray)  and  our  predictions  from  models  Ice  Id  (solid)  and  Ice  2d  (dashed)  at  Krafla,  near 
the  plume  center  (a),  and  at  DSDP  Site  409  on  the  Reykjanes  Ridge  550  km  away  from  the 
plume  center  (b).  This  inversion  method  assumes  fractional  melting  and  includes 
differences  in  partitioning  coefficients  between  the  spinel  and  garnet  stability  fields.  It  also 
assumes  complete  extraction  and  mixing  of  all  melts  generated  in  the  melting  region,  which 
makes  the  estimation  of  maximum  depth  of  melting  sensitive  to  the  low-degree  melt 
compositions  [49].  Another  assumption  is  the  parent  source  composition  (primitive  mantle 
beneath  Krafla  and  a  50-50%  mix  of  primitive  and  depleted  MORE  source  along  the 
Reykjanes  Ridge),  which  is  important  in  estimating  the  maximum  extent  of  melting,  (c) 
Comparison  between  observed  Sr  isotope  concentrations  [41]  along  Iceland  and  the  MAR 
and  weighted  mean  plume  tracer  concentration  P  in  the  accumulated  melts  for  models  Ice 
Id  (solid)  and  Ice  2d  (dashed).  The  peak  in  ^"^Sr/Sesi-  jq  north  of  Iceland  is  due  to  the 
Jan  Mayen  hotspot  [41]  which  we  do  not  model. 
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Fig.  9.  (a)  Lower  diagram  shows  predicted  P-wave  velocity  variations  with  contour 
interval  of  0.5%  for  model  Ice  Id  (broad  plume  source)  caused  by  the  combined  effects  of 
temperature,  melt  depletion,  and  melt  retention.  Top  panel  illustrates  the  predicted  P-wave 
travel-time  delays,  assuming  vertically  passing  rays,  for  along  axial  and  across-axis 
profiles  due  to  successively  added  mantle  effects  as  labeled,  (b)  Same  as  (a)  but  for  model 
Ice  2d  (narrow  plume  source).  The  lowest  velocity  region  occurs  at  depths  50-100  km  due 
to  the  predicted  high  melt  retention. 
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P-wave  Delay  (s) 


Fig.  10.  Characteristic  variables  predicted  by  models  with  melting  normalized  by  those 
predicted  by  models  without  melting.  We  choose  the  mean  value  of  WiQ/U)'^!^  for  each 
experimental  set  and  maximum  value  of  along-axis  variations  for  each  of  the  other 
variables.  Crustal  thickness  anomalies  are  normalized  by  calculations  with  thermal 
buoyancy  and  latent  heat  loss. 


(With  Melting)/(Without  Melting)  (%) 


13 


1— 

o 


Figure  10 


89 


90 


Chapter  4 


Dynamic  Interaction  Between  Mantle  Plumes  and 
Migrating  Midocean  Ridges 


91 


92 


Dynamic  interaction  between  mantle  plumes  and 
migrating  midocean  ridges 


Garrett  Ito 

MIT/WHOI  Joint  Program,  Department  of  Geology  and  Geophysics,  Woods  Hole  Oceanographic 
Institution,  Woods  Hole,  MA,  02543;  508-289-2575  (tel),  508-457-2187  (fax),  gito@magellan.whoi.edu 


Jian  Lin 

Department  of  Geology  and  Geophysics,  Woods  Hole  Oceanographic  Institution,  Woods  Hole,  MA,  02543; 
508-289-2576  (tel),  508-457-2187  (fax),  jlin@whoi.edu 

Carl  W.  Gable 

Earth  and  Environmental  Sciences,  Los  Alamos  National  Laboratories,  Los  Alamos,  NM,  87545;  505-665- 
3533  (tel),  gable@lanl.gov 


Abstract.  We  investigate  the  three-dimensional  interaction  of  mantle  plumes  and 
migrating  midocean  ridges  with  variable  viscosity  numerical  models.  Scaling  laws  derived 
for  stationary  ridges  in  steady  state  with  near-ridge  plumes  are  consistent  with  those 
obtained  from  independent  studies  of  Ribe  [1996].  Our  numerical  results  suggest  that 
along  axis  plume  width  W  and  maximum  distance  of  plume-ridge  interaction  x^ax  scale 
with  where  Q  is  plume  source  volume  flux  and  U  is  ridge  full  spreading  rate. 

Both  W  and  Xmax  increase  with  buoyancy  number  IJt,,  which  reflects  the  strength  of 
gravitational-  versus  plate-driven  spreading,  and  7,  which  is  the  ratio  of  ambient/plume 
viscosity.  In  the  case  of  a  migrating  ridge,  the  distance  of  plume-ridge  interaction  is 
reduced  when  a  ridge  migrates  toward  the  plume  due  to  the  excess  drag  of  the  faster- 
moving  leading  plate,  and  enhanced  when  a  ridge  migrates  away  from  the  plume  due  to  the 
reduced  drag  of  the  slower-moving  trailing  plate.  Thermal  erosion  of  the  lithospheric 
boundary  layer  by  the  plume  further  enhances  W  and  x^ax  but  to  a  degree  that  is  secondary 
to  the  differential  migration  rates  of  the  leading  and  trailing  plates.  These  numerical  models 
are  tested  by  comparing  model  predictions  of  bathymetry  and  gravity  with  observations  of 
the  Galapagos  plume-migrating  ridge  system.  The  amplitudes  and  widths  of  along- 
isochron  bathymetric  and  gravity  anomalies  can  be  explained  with  models  of  a  plume 
source  temperature  anomaly  of  80-120°C,  radius  of  80-100  km,  and  volume  flux  of  4.5  x 
106  km^/m.y.  The  observed  increase  in  anomaly  amplitude  with  isochron  age  is  also 
explained  by  our  models  which  predict  higher  crustal  production  rates  when  the  ridge  was 
closer  to  the  plume  source  several  million  years  ago.  The  same  plume-source  models  also 
predict  crustal  production  rates  of  the  Galapagos  Islands  that  are  consistent  with  those 
estimated  independently  from  the  observed  island  topography.  Predictions  of  the 
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geochemical  signature  of  the  plume  along  the  present-day  ridge  suggest  that  mixing 
between  the  plume  and  ambient  mantle  sources,  as  inferred  from  geochemical  observations, 
is  unlikely  to  occur  in  the  asthenosphere  or  crust.  Instead,  mixing  most  likely  occurs  much 
deeper  in  the  mantle,  possibly  by  entrainment  of  ambient  material  as  the  plume  ascends 
through  the  depleted  portion  of  the  mantle  from  its  deep  source  reservoir. 

INTRODUCTION 

A  wide  range  of  geologic  and  geochemical  observations  provide  strong  evidence  that 
mantle  plumes  feed  material  to  nearby  midocean  ridges  [e.g.  Vogt,  1971;  Schilling,  1973; 
Schilling  et  al,  1976;  Morgan,  1978].  Near-ridge  plumes  are  documented  to  generate 
along-axis  geophysical  anomalies  with  widths  exceeding  2000  km  [Ito  and  Lin,  1995b]  and 
can  induced  geochemical  signatures  for  plume-ridge  separation  distances  of  nearly  1500  km 
[Schilling,  1991].  The  "mantle-plume  source/migrating  ridge  sink"  model  of  Schilling  and 
co-workers  suggests  that  migrating  ridges  are  "fed  and  dynamically  affected  by  a 
preferential  plume  flow  along  a  thermally  induced  channel  at  the  base  of  the  lithosphere" 
[Schilling,  1991].  This  model  suggests  that  a  thermal  channel  is  progressively  carved  into 
the  lithosphere  as  the  ridge  migrates  over  and  away  from  the  impinging  hot  plume 
[Morgan,  1978;  Schilling,  1985;  Schilling  etal,  1985].  All  of  the  13  plume  ridge  systems 
considered  by  Schilling  [1991]  have  ridges  migrating  away  from  their  nearby  plumes  in 
support  of  this  plume  source-migrating  ridge  sink  model. 

Recent  numerical  modeling  and  laboratory  experimental  studies  have  begun  to 
characterize  the  kinematic  and  dynamic  aspects  of  the  interaction  between  mantle  plumes 
and  stationary  midocean  ridges.  For  ridge-centered  plumes,  scaling  laws  for  the 
dependence  of  along-axis  plume  width  W  on  plume  volume  flux  Q  and  ridge  full  spreading 
rate  U  were  first  explored  in  tank  experiments  [Feighner  and  Richards,  1995]  and  further 
developed  in  numerical  studies  [Feighner  et  al.,  1995;  Rite  et  al.,  1995;  Ito  et  al,  1996]. 
The  dynamics  of  off-axis  plumes  were  first  investigated  in  the  laboratory  by  Kincaid  et  al. 
[1995a]  and  in  2-dimensional  (2-D)  numerical  experiments  by  Kincaid  et  al.  [1995b]. 
Finally,  Rite's  [1996]  study  of  off-ridge  plumes  established  scaling  laws  for  the 
dependence  of  VF  on  a  range  of  variables  including  Q,  U,  plume-ridge  distance  Xp,  and 
lithospheric  thickening  with  age. 

While  the  above  studies  established  scaling  parameters  of  plume  ridge  interaction  they 
did  not  investigate  the  effects  of  ridge  migration.  In  the  more  realistic  case  of  a  migrating 
ridge,  not  only  may  thermal  thinning  of  the  lithosphere  be  important  as  proposed  by 


94 


Schilling's  and  co-worker's  plume  source-ridge  sink  hypothesis,  but  also  the  plate  trailing 
the  migrating  ridge  moves  significantly  slower  than  the  plate  leading  the  ridge,  thereby 
inducing  less  drag  on  the  plume  away  from  the  ridge  [Ribe,  1996] . 

We  here  explore  the  dynamics  of  migrating  ridges  and  plumes  with  3-D  numerical 
models  that  include  thermal  diffusion  and  fully  pressure-  and  temperature-dependent 
rheology.  We  will  first  establish  scaling  laws  for  along-axis  plume  width  W  and  maximum 
plume-ridge  interaction  distance  x^ax  for  steady-state  systems  of  stationary  ridges.  These 
results  will  be  compared  with  those  of  the  chemically  buoyant,  constant  viscosity  plume 
models  of  Ribe  [1996]  first  to  verify  the  scaling  parameters  and  then  to  quantify  the 
importance  of  thermal  diffusion  and  variable  plume  viscosity  on  these  scaling  laws.  We 
will  then  quantify  the  effects  of  ridge  migration  on  W  and  Xmax  and  identify  the  physical 
mechanisms  controlling  this  plume  source-migrating  ridge  sink  model.  Finally,  we  will 
test  our  models  by  comparing  model  predictions  with  geophysical  observations  of  the 
Galapagos  plume-migrating  ridge  system,  and  then  discuss  the  implications  for  the 
dimensions,  temperature  anomaly,  volume  flux,  and  geochemical  signature  of  the 
Galapagos  plume. 


GOVERNING  EQUATIONS  AND  NUMERICAL  METHOD 

The  mantle  is  modeled  as  a  viscous  Boussinesq  fluid  of  zero  Reynolds  number  and 
infinite  Prandtl  number.  The  equilibrium  equations  include  conservation  of  mass 


V  •!/  =  0, 

(1) 

momentum 

V  •  T  =  Apg, 

(2) 

and  energy 

K-v^r-H.vr 

dt 

(3) 

(see  Table  1  for  definition  of  variables).  Mantle  density  p 

is  reduced  by  thermal  expansion 

such  that  Ap  =  PoOcAT,  and  the  3-D  stress  tensor  T  depends  on  the  strain  rate  tensor  e 
according  to  T  =  27]  £  -pi.  Viscosity  rj  depends  on  pressure  p  and  real  temperature  Tp 
according  to. 


T]=rio  exp 


E  +  pV 
RTp 


E  +  p„g(0.5D)V 


(4) 
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in  which  Tj^q  is  the  real  temperature  of  the  mantle  at  z  =  0.5D.  Reduced  values  of  R  and  V 
are  used  to  simulate  numerically  the  behavior  of  a  non-Newtonian  rheology  (i.e.  e  ocT  3) 
[Christensen,  1984],  The  ratio  of  ambient/plume  viscosity  7  is  defined  as  rjo/rjp,  where  rip 
is  the  viscosity  of  the  plume  at  z  =  0.5  D. 

Calculations  were  done  using  the  Cartesian  numerical  code  first  written  by  Gable 
[1989]  and  Gable  et  al.  [1991],  and  later  modified  by  Ito  et  al.  [1996]  to  incorporate 
variable  viscosity.  The  numerical  setup  is  illustrated  in  Fig.  1.  Two  spreading  plates  are 
simulated  by  imposing  surface  horizontal  velocities  of  %  =  +UI2  and  %  =  -UI2  on  both 
sides  of  a  model  ridge  axis.  A  plume  is  introduced  by  imposing  a  columnar-shaped 
temperature  anomaly  in  the  lower  portion  of  the  box  at  a  distance  Xp  from  the  ridge  axis. 
The  plume  source  is  hottest  {T  =To  +  ATp)  at  its  center  and  cools  as  a  Gaussian  function  of 
radial  distance  to  Tq  at  its  full  radius.  The  vertical  sides  of  the  box  are  free  of  shear  stress 
and  have  zero  horizontal  temperature  gradient.  Therefore,  the  symmetry  introduced  by  the 
reflecting  side  boundaries  allow  this  half-of  plume-ridge  system  in  numerical  models  to 
simulate  a  full  plume-ridge  system  in  virtual  space.  Temperature  at  the  surface  is 
maintained  at  0°C  which  cools  and  thickens  a  high  viscosity  lithosphere  approximately  with 
the  square  root  of  distance  from  the  ridge  axis.  Temperatures  in  the  lower  portion  of  the 
box  (z  >  0.6D)  are  maintained  at  the  reference  mantle  potential  temperature  Tq  everywhere 
except  inside  the  plume  source.  Correspondingly,  the  energy  equation  is  solved  in  only  the 
upper  portion  of  the  box  (0.6Z)  >  z  >  0).  The  purpose  of  the  lower  volume  of  the  box  is  to 
simulate  an  open  boundary  at  the  base  of  the  upper  volume  where  the  plume-ridge 
interaction  occurs.  To  ensure  numerical  accuracy  of  the  flow  solutions,  we  limit  the 
horizontal  viscosity  variation  to  be  <10^  by  defining  an  upper  viscosity  limit  for  the 
lithosphere  depending  on  the  viscosity  of  the  hot  plume. 

To  track  the  flow  of  the  mantle  plume,  we  introduce  a  passive  tracer  P  in  the  plume 
source  with  value  of  1.0  to  represent  100%  plume  material.  A  finite  difference,  tensor 
diffusion  method  [Gable,  1989;  Travis  et  al.,  1990]  is  used  to  solve  for  advection  of  P, 
from  the  source  and  throughout  the  upper  volume  of  the  box  (z  <  0.6D).  Diffusion  of  P  is 
required  by  our  numerical  method  but  the  rate  of  diffusion  is  reduced  by  a  factor  of  3 
relative  to  the  rate  of  thermal  diffusion.  P  is  also  used  to  determine  along-axis  width  of  the 
plume  by  measuring  the  along-axis  distance  over  which  the  mean  value  of  P  beneath  the 
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ridge, 


(  j  0.6D 


0.6D 


\P{0,y,z)dz 


P  >  0,  is  greater  than  0.1.  Finally,  we  use  P  in  the  steady- 


state  stationary  ridge  case,  to  measure  the  volume  flux  of  plume  material  crossing  the  ridge 
axis  by  integrating  horizontal  velocities  on  the  side  of  the  ridge  opposite  the  plume  where  P 
>  0.4.  The  volume  flux  of  the  plume  source  Q  is  measured  by  integrating  vertical  velocities 
at  z  =  0.6D  over  the  cross-sectional  area  of  the  source  column. 


SCALING  LAWS  FOR  STATIONARY  RIDGES 

Feighner  and  Richards  [1995]  and  Feighner  et  al.  [1995]  demonstrated  that  Wq  = 
is  an  effective  length  scale  for  characterizing  the  horizontal  dimension  of  a  ridge- 
centered,  gravitationally  spreading  plume.  They  also  defined  a  plume  buoyancy  number 
Tib  =  Qollfi  where  a  =  gAp/ASrio,  which  characterizes  the  relative  strength  of  gravitational 
versus  plate-driven  spreading.  Subsequent  analyses  of  Ribe  et  al.  [1995]  derived  a 
characteristic  plume  thickness  scale  So  =  ,  which  determines  Tib  according  to  Tib 

=  (Wo/So)^.  The  effect  of  the  sloping  lithosphere  on  the  interaction  of  off-axis  plumes  was 
characterized  with  the  "upslope  number"  /7„  =  by  Ribe  [1996]. 

The  above  scaling  quantities  were  shown  by  lubrication  theory  models  of  Ribe  [1996] 
to  define  a  full  scaling  law  of  along-axis  width  W  for  steady-state  stationary  ridges, 

w  =  w^Fi(n^)P4(n^,n„)F3(:^,n^,n„).  (5) 

Functions  Fj  and  F4  describe  the  increase  in  steady-state  width  with  increasing  values  of 
lib  and  TIu,  for  ridge-centered  plumes  (xp  =  0);  whereas  function  F3  describes  the  first- 
order  dependence  of  W  on  plume-ridge  separation  distance  Xp  and  the  second-order 
dependence  on  Tib  and  J7„.  We  now  further  investigate  this  scaling  law  with  numerical 
models  that  include  both  thermal  diffusion  and  temperature-dependent  plume  viscosity. 

Ridge-centered  plumes 

The  simplest  case  is  that  of  a  ridge-centered  plume.  In  this  case  Xp  =  0  and  Fj  =  1.0, 
therefore  we  seek  to  define  functions  F]  and  F4.  In  our  numerical  experiments,  we  vary 
full  spreading  rate  U  between  20  and  120  km/m.y.  and  modulate  plume  flux  Q  by  varying 
plume  temperature  anomaly  ATp  between  100  and  200°C  (see  Table  2).  Three  models  of 


97 


plume  viscosity  structure  are  examined.  The  first  is  designed  to  simulate  the  constant 
plume-viscosity  calculations  of  Ribe  [1996].  This  viscosity  structure  omits  the  pressure- 
dependence  of  Eq.  4  and  has  r]  =  rj  ofor  T>To,  thus  plume  viscosity  is  the  same  as  the 
ambient  fluid  (7=  1.0).  To  allow  for  a  thickening  lithospheric  boundary  layer,  we 
incorporate  the  temperature-dependence  in  Eq.  4,  for  T  <  Tq.  The  second  and  third 
viscosity  models  have  the  full  pressure-  and  temperature-dependence  as  defined  by  Eq.  4; 
the  second  has  7  =  2.352  for  ATp  =  100°C,  and  the  third  has  7  =  5.053  for  ATp  =  200°C. 

A  scaling  law  for  normalized  plume  width  WIWo,  which  defines  F;,  is  determined  by 
fitting  WIWo  to  exponential  functions  of  the  quantity  77^7  a  modified  buoyancy  number 
defined  by  the  viscosity  of  the  plume  [Ito  et  al,  1996].  WIWq  is  described  well  by  the 
function 

\ogio(W/Wo)  =  0.32  +  0.01[logio(i7z,7)]  +  0.05  [log  io(77z,7)]2  (6) 

with  a  standard  deviation  misfit  of  8%  of  the  median  value  of  2.25  (Fig.  2).  This  function 
is  consistent  in  general  form  with  Rite's  [1996]  results  of  \o%\q{WIWo)  =  0.217  + 
0.0569[logio(TZ),7)]  +  0.0176[logio(T/^7)]2  -1-  0.0275[logio(i7ft7)]2.  The  relatively  weak 
dependence  of  W  on  logio(TZ^7)  in  our  results  may  reflect  our  source  radius  of  finite 
width,  which  becomes  comparable  to  Wg  at  low  values  of  TZ^7  and  thus  contributes  to 
along-axis  width  in  a  manner  unlike  by  Rite's  [1996]  point  source  plumes.  In  addition, 
we  are  unable  to  identify  a  dependence  on  /7„  which  is  described  by  Rite's  [1996]  function 

=  (1  +  1.77  n„  77^0.33). 

Off-axis  plumes 

To  derive  scaling  laws  for  off-axis  plumes,  we  seek  to  define  the  function  F 3.  Fig.  3 
illustrates  the  shape  of  the  plume  at  different  distances  from  the  ridge  axis.  When  the 
plume  is  ridge-centered,  it  spreads  along  the  ridge-axis,  is  divided  by  the  spreading  plates, 
and  then  spreads  symmetrically  away  from  the  ridge  axis.  When  it  is  off  the  ridge,  the 
plume  spreads  asymmetrically  beneath  the  moving  plate  with  the  upwind  side  tapering 
towards  the  ridge  and  the  downwind  side  widening  away  from  the  ridge  as  it  is  sheared 
away  by  the  moving  plate.  The  ridge  thus  captures  a  narrower  width  of  the  plume  as  Xp  is 
increased.  If  Xp  is  large  enough,  the  ridgeward  flowing  plume  material  stagnates  against 
the  migrating  plate  as  investigated  by  Sleep  [1987]  and  Ribe  and  Christensen  [1994].  It  is 
this  stagnation  distance  that  defines  the  maximum  distance  to  which  plume  material  will 
contact  the  ridge  axis,  Xmax- 
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Function  F3  is  the  dependence  of  W  on  plume-ridge  distance  and  is  equivalent  to 
W/Wixp=0)  (Eq.  5).  Fig.  4a  shows  numerical  results  of  F3  versus  Xp/Wo-  The  best  fitting 
function  is  a  binomial  function  of  the  form 

Fi  =  [  1 .0  -  0.68(xp/WoF2)2]  1/2  (7) 

as  consistent  with  that  of  Ribe  [1996].  As  evident  in  Fig.  4a,  cases  with  7=  1.0  yield  the 
shortest  distances  of  plume-ridge  interaction,  whereas  increasing  values  of  7  result  in 
greater  distances  of  plume-ridge  interaction.  This  second  order  variation  in  plume  width 
reflects  a  stretching  function  F2,  which  depends  primarily  on  7  and  secondarily  on  Tib  with 

a  best  fitting  function 

F2(nbrMnbr)^-^'^r^'^^.  (S) 

As  illustrated  in  Fig.  4b,  F2  collapses  values  of  W/W{xp=0)  on  to  a  single  curve.  Thus  the 
combined  Eqs.  7  and  8  describe  effectively  the  dependence  of  W  on  plume  ridge  distance 
for  steady-state  cases.  The  primary  dependence  of  F2  on  7may  indicate  that  not  only  do 
less  viscous  plumes  spread  stronger  gravitationally,  but  also  they  are  subject  to  less 
shearing  from  the  overlying  migrating  plate.  The  increase  in  Xmax  with  7 as  predicted  here 
is  consistent  with  results  of  2-D  experiment  by  Kincaid  et  al.  [1995b].  F2  derived  from  our 
numerical  models  captures  the  linear  exponential  term  of  Ribe's  [1996]  function, 
logio(J^2)=0.043[logio(nz, )]  +  0.060[logio(nfc  )]2  -  0.0062[logio(nfo  )]3;  however,  our 
results  show  a  weaker  dependence  on  Tib  ■  We  again  do  not  observe  a  strong  dependence 
of  Fj  on  TJu  over  the  range  of  Fr„  examined. 

We  also  investigate  the  percentage  of  the  plume  flux  that  crosses  the  ridge  axis  Qr-  For 
ridge-centered  cases,  half  of  the  plume  material  flows  to  each  side  of  the  ridge  such  that  Qr 
=  0.50.  As  the  plume  moves  away  from  the  ridge  axis,  Qr  decreases  according  to 

Qr  =  0.50  -  0AlixpAVoF2)  (9) 

as  illustrated  in  Fig.  5.  This  function  is  again  similar  to  that  of  Ribe's  [1996],  Qr  =  0.5- 
0.79j!?  +  0.24p^,  where  p  =  {xp/W oF 2)- 

Functions  F3  and  Qr  (Eqs.  7  and  9)  are  zero  when  Xp  =  Xmax,  and  consequently 

Xmax=\-'^^^oF2-  (10) 

For  the  case  of  7  =  1 .0,  our  predicted  values  of  Xmax  l^o  ^re  -50%  greater  than  those 
predicted  by  Ribe  [1996].  Some  of  this  discrepancy  may  be  due  to  differences  in  numerical 
models;  for  example  our  finite  source  radius  versus  Ribe's  [1996]  point  source  as 
mentioned  earlier.  Another  potentially  important  cause  of  this  discrepancy  may  be  thermal 
erosion  of  the  lithosphere.  Fig  6.  illustrates  the  thickness  of  the  lithosphere  that  was 

eroded  by  the  plume  scaled  by  the  modified  characteristic  plume  thickness  = 
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1/4 

.  The  greatest  erosion  occurs  downwind  of  the  plume  where  the  plume  has 

been  in  contact  with  the  lithosphere  the  longest.  The  downwind  slope  of  the  channel  acts  to 
inhibit  ridgeward  spreading  as  noted  by  Kincaid  et  al.  [1995a,b];  however,  our  results 
suggest  that  the  dominant  effect  is  the  slope  of  the  channel  in  the  y-direction  which 
enhances  spreading  toward  the  ridge  by  inhibiting  spreading  in  the  along-axis  direction. 
Ribe  [1996]  predicted  this  effect  to  enhance  Wby  -10%. 

The  similarities  between  the  above  scaling  laws  for  ridge-centered  and  off-axis  plumes 
and  those  of  Ribe  [1996]  indicate  that  the  general  form  of  these  scaling  laws  are  robust  and 
insensitive  to  differences  in  far  field  boundary  conditions.  The  additional  physics  we 
include  are  variable  plume  viscosity  and  thermal  erosion,  both  of  which  enhance  W  at  a 
given  value  of  Xp  as  well  as  the  total  range  over  which  an  off-axis  plume  interacts  with  a 
nearby  stationary  ridge. 

SCALING  LAWS  FOR  MIGRATING  RIDGES 

To  derive  scaling  laws  for  the  case  of  a  migrating  ridge  we  simulate  a  ridge  moving  in 
the  positive  x-direction  at  velocity  VV.  With  respect  to  the  ridge,  both  plates  are  assumed  to 
spread  symmetrically  at  a  rate  of  f//2;  therefore  with  respect  to  the  plume,  Plate  1  (the 
leading  plate  moving  in  the  positive  x-direction)  spreads  with  velocity  +UI2  +  VV  and  Plate 
2  (the  trailing  plate  moving  in  the  negative  x-direction)  spreads  with  velocity  -UH  +  Vr 
(Fig.  7).  These  velocity  conditions  are  incorporated  by  defining  the  appropriate  horizontal 
velocities  of  the  surface  boundary,  whereas  the  motion  of  the  ridge  is  simulated  by 
redefining  the  x-position  of  the  diverging  surface  velocities  at  each  step  during  time 
integration. 

Numerical  experiments  began  with  the  steady-state  configuration  of  a  plume  and 
stationary  ridge,  with  the  plume  beneath  Plate  1  at  x^  >  Xmax-  We  then  allowed  the  ridge  to 
migrate  toward,  over,  and  away  from  the  plume  such  that  the  plume  ends  up  beneath  Plate 
2.  We  use  the  convention,  Xp>0  when  the  plume  is  beneath  Plate  1  and  Xp<0  when  the 
plume  is  beneath  Plate  2.  Three  ridge  migration  velocities  are  tested  for  parameters  of 
experiments  3,  5,  7,  8,  and  12  (Table  2).  In  each  case,  the  maximum  Vr  examined  is  equal 
to  the  half  spreading  rate. 

The  dependence  of  W  on  Xp  is  shown  in  Fig.  8  for  experiment  7.  The  form  of  the 
function  of  W  versus  Xp/Wo  is  the  same  as  that  of  F3  in  Eq.  7,  but  the  curves  are  shifted 
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increasingly  in  the  negative  x^-direction  with  increasing  Vr,  and  the  total  range  over  which 
the  plume  interacts  with  the  ridge  broadens  with  F3  is  thus  modified  to 


F3  = 


lO-OMFeisMbT) 


Xp/W^  +  F^isMbY) 

F2i^Mbr) 


1/2 


S={2VrlU). 


(11) 


Function  F5  describes  the  shift  of  the  curves  in  the  negative  x^-direction  and  is  best  fit  by 
the  function 

F5  =  0.39(ni7)-0l2(2y^f/)  (12) 

(Fig.  9a).  Function  F^  controls  the  increase  in  total  range  of  plume-ridge  interaction  and  is 
best  fit  by  the  function 

Fg  =  1.0  -  0.17(i7/,7)-O  l2(2y^[/)l/2  (13) 

(Fig.  9b). 

Function  F5  reflects  largely  the  differential  shearing  of  the  plume  by  the  asymmetrically 
moving  plates  and  has  the  largest  effect  on  the  range  of  plume  ridge  interaction.  When  Xp  > 
0,  the  plume's  upwind  stagnation  point  defines  the  maximum  distance  to  which  the  plume 
interacts  with  the  ridge.  The  faster  moving  Plate  1  induces  more  drag  on  the  plume  away 
from  the  ridge  (Fig  7a)  therefore  pushing  the  stagnation  point  closer  to  the  plume  source 
and  reducing  Xmax  relative  to  the  case  in  which  Vr  =  0.  When  Xp  <  0,  the  plume  separates 
from  the  ridge  when  the  ridgeward  spreading  velocity  of  the  plume  drops  below  the 
migration  rate  of  the  ridge.  The  slower  moving  Plate  2  induces  less  shear  away  from  the 
ridge  (Fig.  7b),  consequently,  the  plume  is  able  to  keep  up  with  the  migrating  ridge  over  a 
greater  distance.  F5  reduces  Xmax  for  x  >  0  and  increases  Xmax  for  x  <  0  by  as  much  as 
35%  of  Xmax  for  a  stationary  ridge.  The  degree  to  which  the  differential  motion  of  the  two 
plates  is  able  to  alter  the  shape  of  a  plume  diminishes  for  strong  plumes.  This  is  reflected 
in  the  inverse  relationship  between  F5  and  77^7- 

Function  F^  most  likely  reflects  the  effects  of  lithospheric  erosion,  which  increases  the 
total  range  over  which  the  plume  interacts  with  the  ridge.  Thermal  erosion  has  the 
strongest  effect  on  the  system  after  the  ridge  has  migrated  over  and  away  from  the  plume 
(Fig  10)  as  hypothesized  by  Schilling  [1991].  In  the  case  that  Vr  is  small  (Fig.  10a),  the 
velocity  of  Plate  2  is  largest  thus  allowing  the  plume  to  erode  a  greater  area  of  the  plate 
downwind  of  plume  source.  Consequently,  the  plume  spreads  more  easily  away  from  the 
ridge  through  the  downwind  eroded  channel.  On  the  other  hand,  when  Vr  is  larger  (Fig. 
10b  and  10c),  Plate  2  moves  slower,  and  the  downwind  channel  is  more  confined  to  the 
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plume  source,  thereby  enhancing  spreading  toward  the  ridge.  The  dependence  of  on 
2VrlU  weakens,  however,  with  increasing  TlbY  because  the  shape  of  the  overlying 
lithosphere  has  less  effect  on  gravitational  spreading  of  the  plume  at  higher  values  of  77^7 . 
For  Vr  =  UI2,  Fb  increases  the  total  range  of  plume  ridge  communication  by  an  average  of 
11%. 

Thus,  lithospheric  erosion  has  only  second-order  importance  in  influencing  the  flow  of 
near-ridge  plumes — a  result  that  differs  from  that  envisioned  in  Schilling  and  co-worker's 
plume  source-migrating  ridge  sink  model.  The  effects  of  the  lithosphere  are  likely  to  be 
weakest  when  the  characteristic  plume  thickness  is  large  relative  to  the  thickness  of 

the  lithospheric  boundary  layer.  Indeed,  the  values  of  SoY'^'^  examined  here  are  100-150 
km — 3-5  times  the  thickness  of  the  lithosphere  overlying  the  plumes.  The  regime  in  which 
SoY'^^'^  is  comparable  to  the  thickness  of  the  lithosphere  would  allow  the  lithosphere  to 
influence  more  strongly  the  ability  of  the  plume  to  spread  to  the  ridge.  This  low-^o/'^^^ 
regime  would  require  significantly  hotter  plumes  to  reduce  rjp,  as  well  as  significantly 
narrower  sources  radii  to  limit  Q.  Such  conditions,  however,  may  be  unusual  in  the  Earth 
given  that  the  50-km  radius  and  100-200°C  temperature  anomalies  examined  here  are 
reasonable  properties  of  Earth  plume  examples  [e.g.  Ito  and  Lin,  1995b;  Schilling,  1991; 
Wolfe  etal,  1996] 

Our  complete  scaling  law  for  plume  width  and  migrating  ridge  is  thus 


W=WoF](nbY) 


l.0-0.6SF^(s,nbY) 


Xp/W^  +  F5(s,nbY) 

Fli^MbY) 


^2 


il/2 


(14) 


and  our  complete  our  scaling  law  for  the  maximum  plume-ridge  interaction  distance  is 

X,nax  =  i±l.2\F2F6-F2.F5)Wo.  (15) 

Ridge  migration  has  first  order  effeets  on  the  dynamics  of  plume-ridge  interaetion.  On 
average,  ridges  migrating  toward  plumes  at  rates  comparable  to  their  half  spreading  rates, 
can  sample  plume  material  over  distances  -24%  less  than  stationary  ridges.  Ridges 
migrating  away  from  plumes  at  rates  comparable  to  their  half  spreading,  however,  are  able 
to  sample  plume  material  to  plume-ridge  distances  -36%  greater  than  stationary  ridges  and 
almost  twice  as  far  as  ridges  migrating  toward  plumes. 


THE  GALAPAGOS  PLUME-MIGRATING  RIDGE  SYSTEM 
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We  compare  model  predictions  with  observations  at  the  Galapagos  plume  and 
spreading  center,  a  classic  and  relatively  well  studied  example  of  an  off-axis  plume- 
migrating  ridge  system  (Fig.  11).  The  Galapagos  spreading  center  separates  the  Cocos 
plate  to  the  north  and  the  Nazca  plate  to  the  south  with  a  full  spreading  rate  of  ~55  km/m.y. 
at  9rW  [DeMets  et  al.,  1994].  The  ridge  is  currently  migrating  northward  with  respect  to 
the  hotspot  at  a  rate  of  27  km/m.y.  [Gripp  and  Gordon,  1990],  which  is  ~1  km/m.y.  less 
than  the  half-spreading  rate,  ho  and  Lin  [1995a]  documented  that  the  total  amplitude  of 
bathymetric  and  mantle-Bouguer  gravity  anomalies  (MBA)  along  Cocos  Plate  isochrons 
increase  with  isochron  age,  and  suggested  that  this  behavior  reflects  increased  crustal 
production  in  the  past  when  the  plume  was  closer  to  the  ridge.  Our  purpose  here  is  to 
compare  predicted  and  observed  profiles  of  bathymetry  and  MBA  in  order  to  assess  the 
degree  to  which  the  models  can  explain  the  observations  and  to  place  theoretical  constraints 
on  the  dimensions,  temperature  anomaly,  and  flux  of  the  Galapagos  plume. 

Calculations  of  crustal  thickness,  bathymetric,  and  gravity  anomalies 

Because  70-75%  of  the  along-isochron  bathymetric  and  gravity  variations  most  likely 
arise  from  plume  induced  thickening  of  the  igneous  crust  [Ito  and  Lin,  1995a],  crustal 
thickness  calculations  are  a  crucial  link  between  our  fluid  dynamic  models  and  surface 
observations.  To  predict  crustal  thickness  along  a  model  ridge  axis,  we  incorporate  the 
solidus  and  liquidus  functions  of  McKenzie  and  Bickle  [1988],  as  well  as  their  functional 
dependence  of  melt  fraction  M  on  homologous  temperature  for  adiabatic  batch  melting. 
Assuming  melt  generated  in  the  mantle  accretes  perpendicularly  to  the  ridge  axis,  crustal 
thickness  along  the  ridge  is  calculated  according  to 

Cr{y)  =  —  —  |[M(x,y,z)dxdz  (16) 

UyPmr 

where  M  =  .  This  method  generates  a  normal  ridge  crustal  thickness  of  6.5  km 

dt 

with  the  assumed  ambient  mantle  temperature  Tq  of  1300  °C.  Because  the  Galapagos 
plume  enhances  crustal  production  at  the  ridge  as  well  as  generates  Galapagos  Islands,  an 
important  source  of  uncertainty  is  how  melt  produced  by  the  plume  is  partitioned  between 
the  ridge  and  hotspot  islands.  We  do  not  attempt  to  model  melt  migration  from  the  mantle 
to  the  ridge  and  islands,  but  instead,  we  assume  that  all  melt  generated  closest  to  the  ridge 
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axis  accretes  at  the  Galapagos  Spreading  Center  and  all  melt  generated  closest  to  the  plume 
source  accretes  at  the  Galapagos  Islands.  Crustal  thickness  along  the  ridge  and  at  the 
hotspot  therefore  change  with  time  as  a  direct  result  of  the  position  of  the  ridge  and  plume 
source.  For  this  reason  we  compare  not  only  predicted  and  observed  anomalies  at  the 
Galapagos  Spreading  Center  but  also  crustal  production  rates  of  the  Galapagos 
Archipelago. 

When  considering  melting,  it  is  important  also  to  account  for  its  effects  on  the  mantle 
[ho  et  al.  1996].  Melting  reduces  mantle  temperature  due  to  latent  heat  loss,  which 
increases  both  mantle  density  and  viscosity;  but  at  the  same  time,  melting  reduces  mantle 
density  by  preferential  extraction  of  iron  with  respect  to  magnesium  [Oxburgh  and 
Parmentier,  1977].  Latent  heat  loss  is  incorporated  by  introducing  a  source  term  - 
{TAS/cp)M  in  the  energy  equation  (Eq.  3).  The  compositional  effect  on  mantle  density  is 
incorporated  by  the  equation 

Ap  ^  PoiaT  +  pX),  (17) 

where  X  is  the  extent  of  melt  depletion  and  =  0.24  is  a  coefficient  of  depletion  density 
reduction  [Oxburgh  and  Parmentier,  1977].  The  equilibrium  equation  for  the  depletion 
field  is 

=  -u.XX  +  M,  (18) 

at 

in  which  the  advection  term  is  solved  using  the  same  tensor  diffusion  method  as  that  used 
to  solve  for  temperature  field,  and  the  source  term  M  is  solved  as  describe  above.  The 
above  melting  effects  do  not  modify  significantly  the  broad  scale  flow  of  the  plume  [Ito  et 
al.  1996];  however,  they  contribute  substantially  to  the  mantle  contributions  to  bathymetric 
and  gravity  anomalies. 

To  calculate  isostatic  topography  of  the  seafloor,  we  consider  contributions  from  both 
the  crust  and  mantle.  In  calculating  crustal  topography,  we  assume  Airy-type 
compensation  of  the  crust  assuming  a  normal  crustal  density  of  2700  kg/m^  that  increases 
linearly  along  axis  to  Pmax  within  500  km  of  the  point  closest  to  the  hotspot  (~91°W). 
Values  of  pmax  considered  are  2900  and  3000  kg/m^.  In  calculating  topography  due  to  the 
mantle,  we  assume  Pratt-type  compensation  with  a  compensation  depth  of  200  km  and 
include  both  thermal  and  compositional  density  effects  as  defined  in  Eq.  17. 

Mantle-Bouguer  gravity  anomaly  (MBA)  is  the  free-air  gravity  minus  the  attraction  due 
to  topography  of  the  seafloor  and  crust-mantle  interface  assuming  a  reference  crust  of 
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uniform  density  (e.g.,  2700  kg/m^)  and  thickness  (e.g.,  6.5  km)  [e.g.  Kuo  and  Forsyth, 
1988;  Lin  et  al,  1990],  MBA  therefore  reflects  crustal  thickness  stmcture  that  differs  from 
this  reference  crust  as  well  as  variations  in  mantle  density.  To  calculate  MBA,  we  again 
include  the  contribution  of  along-axis  crustal  thickness  variations,  and  thermal  and 
compositional  mantle  density  variations  {Ito  etal.  1996]. 

We  investigate  two  radii  and  temperature  anomalies  for  the  Galapagos  plume  source: 
one  has  a  radius  of  100  km  and  temperature  anomaly  ATp  of  80°C,  and  the  other  has  a 
radius  of  80  km  and  ATp  of  120°C.  Both  plume  source  models  predict  comparable  volume 
fluxes  of  4.5x1  O^km^/m.y. — a  value  slightly  greater,  but  comparable  to  the  prediction  of 
2.6-3.6  km^/m.y.  by  Schilling  [1991]  and  the  2.2x10^  km^/m.y.  lower-bound  prediction 
of  Ito  and  Lin  [1995b].  Values  for  To  of  1300  °C  and  of  3  x  10^9  Pa  s  yield  a  Rayleigh 
number  of  3.05  x  10^.  We  began  model  calculations  began  with  a  steady-state  condition  of 
the  plume  beneath  Plate  1  (Cocos  Plate).  We  then  activated  ridge  migration  and  tracked 
crustal  production  and  mantle  evolution  as  the  ridge  migrated  over  the  plume  source. 
Calculations  finished  with  Plate  2  (the  Nazca  Plate)  over  the  plume  source  and  when  the 
ridge  was  200  km  from  the  plume  source.  This  distance  is  the  average  distance  between 
Fernadina  Island  and  the  two  ridge  segments  east  and  west  of  the  transform  fault  at  91°W 
(Fig  11). 

Predicted  and  observed  bathymetric  and  gravity  anomalies 

Model  predictions  of  bathymetry  and  MBA  are  compared  with  five  along-isochron 
profiles  investigated  by  Ito  and  Lin  [1995a]:  the  present-day  ridge  axis  and  isochrons  at 
crustal  ages  of  2.6,  3.6,  6.0,  6.6,  and  7.7  m.y.  Hey  [1977]  suggested  that  the  Galapagos 
Spreading  Center  was  centered  over  the  plume  ~10  Ma.  Therefore,  we  associated  model 
crustal  profiles  with  isochrons  by  taking  the  crustal  profiles  generated  with  Xp  values  that 
were  the  same  fractions  of  200  km  as  each  isochron  was  of  10  m.y.  For  example  the  3.6 
Ma  isochron  was  assumed  to  have  formed  when  the  plume  was  36%  closer  to  the  ridge 
than  it  is  today,  i.e.,  =  128  km.  Mantle  bathymetric  and  gravity  profiles  were  extracted 

at  x-positions  on  Plate  1  (Cocos  Plate)  corresponding  directly  to  the  isochron  ages.  The 
total  predicted  bathymetric  and  MBA  anomaly  profiles  are  the  sum  of  the  crustal  and  mantle 
contributions. 

The  comparisons  between  model  and  observed  profiles  in  bathymetry  are  shown  in 
Fig.  12  for  pmax  =  2900  kg/m^.  Along  the  present  day  ridge  axis  and  isochrons  younger 
than  6  m.y.,  both  models  predict  reasonably  well  the  amplitudes  and  wavelengths  of  the 
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observations.  Along  isochrons  older  than  6  m.y.,  the  cooler  plume  source  of  model  1  also 
yields  predictions  consistent  with  the  observations  but  the  hotter  plume  source  of  model  2 
over-predicts  the  bathymetric  anomalies. 

The  similar  anomaly  amplitudes  predicted  by  the  two  models  at  the  youngest  isochrons 
(i.e.,  200  km  >Xp>  128  km)  reflects  the  similarity  between  the  amount  of  melt  partitioned 
to  the  axial  crust  despite  the  differences  in  plume  source  properties.  Although  the  cooler 
plume  source  of  model  1  is  predicted  to  generate  less  total  melt  than  the  hotter  source  of 
model  2,  the  greater  radius  of  the  model  1  plume  source  causes  more  melting  to  occur  near 
the  ridge  axis,  thus,  a  larger  percentage  of  the  total  melt  liberated  is  partitioned  to  the  ridge. 
On  the  other  hand,  the  narrower  source  of  model  2  predicts  melting  to  be  more  localized  to 
the  center  of  the  plume  source,  therefore,  a  smaller  percentage  of  the  total  melt  generated  is 
partitioned  to  the  ridge  axis.  This  trade-off  between  source  radius  and  temperature  explains 
why  at  the  younger  isochrons,  the  two  different  plume  sources  yield  similar  crustal 
thicknesses  at  the  ridge  axis.  Along  the  oldest  three  isochrons  (i.e.,  128  km  >  Xp>  40 
km),  however,  the  differences  between  the  bathymetric  predictions  of  models  1  and  2  are 
greatest  because  the  amount  of  melt  partitioned  to  the  ridge  crust  reflects  a  larger  percentage 
of  the  total  melt  produced.  Consequently,  the  hotter  source  model  (model  2)  over  predicts 
the  crustal  thickness  at  the  ridge  axis. 

While  the  difference  between  the  two  source  temperature  anomalies  is  but  slight,  the 
differences  between  predicted  crustal  thickness  anomalies  at  the  three  oldest  isochrons  is 
substantial:  model  2  predicts  crustal  thickness  anomalies  of  1 1-15  km,  about  twice  as  large 
as  the  6-8  km-anomalies  predicted  by  model  1.  These  large  differences  in  predicted  crustal 
thickness  anomalies  reflects  the  high  sensitivity  of  upwelling  thus  melting  rate  to  plume 
temperature  anomaly  due  to  the  combined  effects  of  reduced  viscosity  and  enhanced 
thermal  and  depletion  buoyancy.  Directly  above  the  plume  source,  model  2  predicts  a 
maximum  upwelling  rate  250  km/m.y.  This  rate  is  nearly  twice  as  high  as  that  predicted  by 
model  1  of  140  km/m.y.,  thus  explaining  the  factor  of  two  difference  between  the  model  2 
and  model  1  crustal  thickness  anomalies  at  the  three  oldest  isochrons. 

Directly  beneath  the  plume-affected  portion  of  the  present-day  ridge  axis  (i.e.  directly 
north  of  the  plume),  however,  both  models  1  and  2  predict  a  -30%  reduction  of 
upwelling/melting  rate  relative  to  that  beneath  the  unaffected  portions  of  the  ridge;  the 
reason  being,  is  that  upwelling  that  normally  accommodates  plate  spreading  is  replaced  by 
lateral  flow  supplied  by  the  plume.  A  significant  proportion  the  melting  that  contributes  to 
ridge  axis  crust  at  the  present-day,  is  therefore  predicted  to  occur  near  the  midpoint  between 
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the  plume  source  and  the  Galapagos  Spreading  Center.  The  maximum  extent  of  melting 
predicted  at  present-day  is  23%  for  model  1  and  26%  for  model  2. 

Because  crustal  thickness  at  the  ridge  is  predicted  to  increase  with  isochron  age,  the 
contribution  to  bathymetry  of  the  crust  relative  to  that  of  the  mantle  also  is  predicted  to 
increase.  Model  1,  for  example,  predicts  a  crustal  uplift  at  the  present-day  ridge  axis  of  0.6 
km,  which  is  60%  of  the  predicted  total  anomaly  of  1.0  km.  In  contrast,  along  the  7.7- 
m.y.  isochron,  model  1  predicts  a  crustal  uplift  of  1 .2  km,  which  is  80%  of  the  predicted 
total  bathymetric  anomaly  of  -1.5  km.  Likewise,  model  2  predicts  a  crustal  uplift  along  the 
present-day  ridge  of  0.6  km,  which  is  50%  of  the  total  anomaly  of  1.2  km,  and  a  crustal 
uplift  along  the  7.7-m.y.  isochron  of  2.4  km,  which  is  80%  of  the  total  predicted  anomaly 
of  -3.1  km.  These  predictions  are  consistent  with  the  gravity  and  bathymetry  analyses  of 
Ito  and  Lin  [1995a]  which  suggested  that  the  depth  of  compensation  shallows  with 
increasing  isochron  age. 

Comparisons  between  predicted  and  observe  MBA  profiles  are  shown  in  Fig  13. 
Similar  to  the  results  of  the  bathymetry  comparisons,  both  models  1  and  2  yield  MBA 
amplitudes  and  widths  consistent  with  the  observations  for  the  three  youngest  isochrons, 
but  model  2  over-predict  the  amplitudes  of  the  MBA  at  the  three  oldest  isochrons.  In  model 
1,  the  crustal  component  of  MBA  is  predicted  to  be  65%  of  the  total  predicted  anomaly  of  - 
80  mGal  at  the  present-day  ridge  axis  and  -82%  of  the  -140-mGal  anomaly  predicted  along 
the  7.7-m.y.  isochron.  Lateral  density  variations  in  the  mantle  supply  the  remaining 
proportions  of  the  anomalies.  In  model  2,  the  crustal  contribution  to  MBA  is  predicted  to 
be  55%  of  the  total  predicted  anomaly  of  -94  mGal  along  present-day  ridge  axis  and  80% 
of  the  total  predicted  anomaly  of  -261  mGal  along  the  7.7-m.y.  isochron.  Thus,  for  both 
along-isochron  MBA  and  bathymetric  anomalies,  the  relative  contribution  of  the  crust  is 
predicted  to  be  50-80%  of  the  total  anomalies — a  range  slightly  greater  than  the  estimates  of 
Ito  and  Lin  [1995a]  who  used  a  passive  mantle  upwelling  model. 

The  predicted  and  observed  amplitudes  of  along-isochron  MBA  and  bathymetric 
anomalies  for  both  pmax  =  2900  and  3000  kg/m^  are  plotted  versus  isochron  age  in  Fig  14 
(a)  and  (b).  For  the  youngest  isochrons,  the  observed  amplitudes  appear  to  be  matched 
best  by  predictions  of  the  hotter  source  of  model  2  with  the  upper-bound  pmax  of  3000 
kg/m3.  For  the  oldest  isochrons,  the  observed  anomaly  amplitudes  are  best  matched  by  the 
cooler  source  of  model  1,  but  again  with  pmax  of  3000  kg/m^.  Model  2  yields  upper- 
bound  predictions  for  the  oldest  isochrons.  In  general,  the  observed  anomaly  amplitudes 
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are  bracketed  by  the  predictions  resulting  from  the  range  of  p^ax  as  well  as  source  radii 
and  temperatures  considered. 

Galapagos  Archipelago  crustal  volume  flux 

As  discussed  above,  because  of  the  large  uncertainty  in  evaluating  how  melt  is 
partitioned  between  the  ridge  crust  and  hotspot  islands  we  must  consider  also  the  crustal 
production  rate  at  the  Galapagos  Archipelago.  We  first  estimate  the  total  volume  of  the 
Galapagos  Archipelago  by  assuming  the  observed  bathymetry  is  supported  by  Airy 
isostasy  of  the  crust.  We  consider  the  bathymetry  in  the  white  box  in  Fig.  11,  the 
longitudinal  extent  of  which  corresponds  to  ~10  m.y.  of  island  accretion  [Sinton  et  al, 
1996].  Lithospheric  flexure  [Feighner  and  Richards,  1994]  is  neglected  here  because 
flexure  acts  to  only  smooth  topography  of  the  crust-mantle  interface  but  does  not  affect  the 
total  volume  of  the  compensating  crustal  root.  To  correct  for  the  long  wavelength  swell 
topography,  which  is  unlikely  to  reflect  island  volcanism,  we  subtract  a  reference  depth 
plane  with  the  box's  average  bathymetric  slope  in  both  longitudinal  and  latitudinal 
directions.  From  this  residual  bathymetric  map,  we  calculate  the  isostatic  thickness  of  the 
Galapagos  Archipelago  and  then  integrate  along  latitudinal  profiles  to  derive  excess  crustal 
volume  as  a  function  of  longitude  across  the  box  (the  mean  thickness  of  a  normal  oceanic 
crust  of  6.5  km  is  excluded).  Each  longitude  is  then  assigned  an  age  assuming  a  constant 
eastward  migration  rate  of  the  Nazca  Plate  relative  to  the  plume.  Finally,  we  estimate 
crustal  volume  flux  as  a  function  of  age  by  dividing  the  estimated  volumes  by  the  time 
spans  represented  by  their  spacing  in  longitude.  To  be  consistent  with  the  assumed  values 
of  Pmax  along  the  ridge  axis,  we  consider  island  crustal  densities  of  2900  and  3000  kg/m^. 

Fig.  14c  shows  the  estimated  island  fluxes  through  time  which  yield  averages  of  1.2  x 
10^  and  1.6  x  10^  km^/m.y.  over  the  past  7.7  m.y.  for  crustal  densities  of  2900  and  3000 
kg/m^  respectively.  Similar  to  the  comparisons  of  the  isochron  anomalies,  the  hotter  plume 
source  in  model  2  predicts  an  island  crustal  flux  most  consistent  with  the  calculated  fluxes 
over  the  most  recent  4  m.y.  and  an  upper-bound  for  the  island  flux  at  times  >  4  Ma.  The 
cooler  plume  source  of  model  1,  on  the  other  hand,  predicts  lower-bound  island  fluxes 
over  the  most  recent  4  Ma  and  fluxes  that  are  more  consistent  with  the  estimated  fluxes  at 
times  >  4  Ma.  In  general,  the  range  of  source  temperatures  and  radii  considered  by  our  two 
models  yield  island  fluxes  consistent  with  those  estimated  from  the  bathymetry  of  the 
Galapagos  Archipelago  . 
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It  is  possible  that  the  Galapagos  plume  source  may  have  changed  through  time  in 
temperature  anomaly,  radius,  or  both  as  hinted  by  the  closer  match  of  the  hotter  source 
model  to  observations  associated  with  crustal  ages  <  4  Ma  and  closer  match  of  the  cooler 
source  model  to  observations  associated  with  crustal  ages  >  4  Ma.  However,  given  the 
range  of  uncertanties  in  our  models  it  is  impossible  to  resolve  such  changes  in  source 
properties.  The  conclusions  we  make  are  that  our  numerical  plume-ridge  models  are 
capable  of  explaining  the  first  order  variations  in  ridge-axis  anomalies  and  island  flux 
estimates  at  present-day,  as  well  as  explaining  the  apparent  evolution  over  the  past  ~8  m.y. 

A  potentially  important  test  of  our  models  would  be  a  mantle  teleseismic  study  of  the 
Galapagos  plume-ridge  system,  which  would  test  directly  our  predictions  of  source 
dimension  and  temperature  anomaly.  Beneath  the  Galapagos  Archipelago,  we  predict  a  P- 
wave  velocity  reduction  of  0.5-0.7%  due  to  the  excess  temperature  of  the  plume  and  up  to 
2%  in  the  melting  region  if  there  is  up  to  3%  of  melt  present  in  the  mantle  [Ito  et  al,  1996]. 
This  prediction  is  based  on  a  6.25x10-3%  reduction  of  P-wave  velocity  for  each  1°C 
temperature  anomaly  and  a  1.25%  decrease  in  velocity  for  each  1%  porosity  of  melt  in  the 
mantle  [Humphreys  and  Dueker,  1994].  Such  velocity  anomalies  are  predicted  to  result  in 
a  0.3-0.4  s  delay  over  the  center  of  the  hotspot  for  P-waves  passing  vertically  through  the 
upper  400  km  of  mantle  we  have  modeled.  Along  the  Galapagos  Spreading  Center, 
however,  we  predict  mantle  P-wave  velocities  to  actually  increase  by  up  to  0.5%  in  the 
melting  zone  relative  to  normal  ridge  mantle.  The  reason  for  this  velocity  increase  is  that 
the  plume  material  feeding  the  ridge  has  already  experienced  melting  at  the  hotspot; 
consequently,  the  velocity  enhancing  effects  of  melt  depletion  (0.1%  velocity  increase  for 
each  1%  degree  of  depletion  [Humphreys  and  Dueker,  1994])  dominate  over  the  velocity 
reducing  effects  of  temperature  and  melt  retention  directly  beneath  the  ridge.  Another 
valuable  study  would  be  to  obtain  seismic  constraints  on  crustal  thickness  variations  along 
the  ridge  axis  and  along  the  seafloor  isochrons.  This  study  would  test  directly  our 
predictions  of  along-isochron  crustal  thickness  variations  and  place  hard  constraints  for 
geodynamic  models  such  as  these. 

Geochemical  implications 

Much  of  the  original  observations  that  led  to  the  concept  that  plumes  feed  nearby  ridges 
comes  from  systematic  variations  in  basalt  chemistry.  Schilling  and  co-workers  noted  that 
Galapagos  ridge  axis  basalts  erupted  nearest  the  Galapagos  hotspot  have  compositional 
affinities  to  ocean  island  basalts  (OIB) — being  enriched  relative  to  midocean  ridge  basalts 
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(MORE)  in  radiogenic  isotopes  [Verma  and  Schilling,  1982;  Verma  et  ah,  1983]  and 
incompatible  rare-earth  and  major  elements  [Schilling  et  al,  1976;  Schilling  et  al,  1982], 
They  showed  that  the  OIB  signatures  decrease  along  the  ridge  axis  with  increasing  distance 
from  the  hotspot.  An  example  of  this  behavior  is  revealed  in  La/Sm  ratios  as  shown  in  Fig. 
15.  Such  a  systematic  decrease  in  the  OIB  signature  is  interpreted  to  reflect  mixing 
between  the  OIB  plume  source  with  the  MORE  upper  mantle  source  material. 

To  investigate  the  processes  of  plume  and  ambient  mantle  mixing,  we  calculate  the 
amount  of  plume  tracer  P  composing  the  model  crust  along  the  ridge.  After  Ito  et  al. 
[1996],  the  average  plume  tracer  concentration  in  accumulated  melts  as  a  function  of  along- 
axis  coordinate  is 


_  f  R(x:,  y,  z)M{x,  y,  z)dxdz 

P{y)  =  - — n - •  (19) 

J  M{x,y,z)dxdz 

By  the  definition,  P  =  1.0  indicates  that  all  melts  generated  in  a  plane  perpendicular  to  that 
point  of  the  ridge  is  entirely  plume-source  derived.  Likewise,  P  =  0.0  indicates  that  none 
of  the  melts  are  plume  derived,  and  0.0  <  P  <1.0  indicates  some  of  the  melts  are  derived 
from  the  plume  and  some  are  derived  from  the  ambient  mantle  material. 

As  shown  in  Fig.  15,  both  models  1  and  2  predict  geochemical  plume  widths  consistent 
with  the  ~1000-km  width  inferred  from  the  La/Sm  anomaly.  The  largest  difference 
between  predicted  and  observed  profiles  is  that  the  predicted  profiles  indicate  very  little 
mixing  between  the  plume-derived  and  ambient  mantle-derived  melts  over  most  of  the 
plume-affected  portions  of  the  ridge  axis.  Only  at  the  outermost  -200  km  within  the  edges 
of  the  plume  is  there  evidence  for  plume-ambient  source  mixing  in  our  models.  Similar  to 
Ito  et  al.'s  [1996]  conclusions  for  Iceland,  we  suggest  that  mixing  does  not  occur  in  the 
shallow  mantle  or  in  the  cmst  but  most  likely  deeper  in  the  mantle  than  we  have  considered 
in  our  models.  Such  a  deep  mixing  process  may  be  entrainment  of  the  ambient  mantle 
material  by  the  plume  as  it  ascends  through  the  isotopically  depleted  region  of  the  mantle 
[e.g.,  Geist  et  al.,  1988;  Graham  et  al.,  1993]. 

DISCUSSION 

The  above  comparisons  of  predictions  and  observations  at  the  Galapagos  system  as 
well  as  the  scaling  laws  for  W  and  Xmax  assume  that  the  along-axis  bathymetric,  MBA,  and 
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geochemical  anomalies  reflect  directly  the  width  of  the  plume  in  the  mantle.  This  is  likely 
the  case  if  melt  migration  along  the  ridge  axis  is  small  or  non-existent.  If,  however,  along- 
axis  melt  migration  is  significant  as  suggested  by  Ito  et  al.  [1996]  for  the  Iceland-Mid- 
Atlantic  Ridge  system,  then  plume  widths  defined  from  geophysical  or  geochemical 
observations — which  reflect  largely  the  properties  of  the  accumulated  crust — are  likely  to 
be  broader  than  the  width  of  the  plume  in  the  mantle.  If  this  is  the  case,  then  the  same 
values  of  W  and  possibly  x^ax  examined  here  may  require  smaller  values  of  Q  than 
suggested  by  our  scaling  laws.  The  implication  for  the  Galapagos  system  is  that  the 
Galapagos  plume  source  may  be  hotter  and  narrower  than  what  our  models  imply. 

Additional  complexities  that  may  affect  the  systematics  of  along-axis  plume  width  and 
Xmax  at  plume-migrating  ridge  systems  are  ridge  jumps  and  asymmetric  plate  spreading. 
Episodes  in  which  the  ridge  jumps  toward  the  neighboring  plume  has  been  documented  for 
the  Galapagos  system  [Wilson  and  Hey,  1995]  as  well  as  other  systems  in  the  southern 
ocean  [Small,  1995].  Such  episodes  may  result  directly  from  plume-ridge  interaction  as  the 
plume  weakens  the  overlying  plate  near  the  ridge  [Small,  1995].  Asymmetrically  spreading 
ridges,  which  may  also  result  directly  from  plume  weakening  of  the  lithosphere,  are  also 
common  to  plume-ridge  systems  [Small,  1995].  Factors  such  as  these  that  affect  the 
relative  motion  of  the  ridge  are  likely  to  have  little  affect  on  x^ax  when  the  ridge  migrates 
toward  the  hotspot  because  in  this  case  x^ax  is  controlled  by  the  stagnation  point  of  the 
plume  rather  than  motion  of  the  ridge.  On  the  other  hand,  ridge  jumps  and  asymmetric 
spreading  may  increase  x^ax  significantly  when  a  ridge  migrates  away  from  the  hotspot 
because  in  this  case  x^ax  is  determined  by  the  point  at  which  the  migrating  ridge  outruns 
the  ridgeward  spreading  plume. 

Regardless  of  how  plume  material  is  sampled  by  midocean  ridges,  our  numerically 
derived  scaling  laws  suggest  that  plumes  affect  broad  regions  of  oceanic  plates.  In  general, 
Eq.  14  and  15  suggest  that  the  maximum  along-axis  width  of  a  plume  is  125-200%  as 
broad  as  the  maximum  plume-ridge  interaction  distance.  The  major  implication  is  that — as 
in  the  Atlantic  and  southern  oceans  with  documented  plume  signatures  at  ridges  located  as 
far  as  1400  km  away — individual  plumes  may  spread  over  distances  of  up  to  2500  km 
perpendicular  to  the  direction  of  plate  motion.  Such  ridge-perpendicular  spreading  may 
generate  broad  bands  of  plume-affected  lithosphere,  which  may  alter  otherwise  normal 
lithosphere  and  contribute  to  characteristic  properties  of  "tectonic  corridors"  such  as  those 
identified  by  Kane  and  Hayes  [1992]  and  Hayes  and  Kane  [1994].  Among  the  most 
prominent  examples  of  plume  affected  lithosphere  are  the  broad  regions  of  anomalously 


111 


shallow  seafloor  associated  with  the  Galapagos  system  as  discussed  in  this  study,  the 
Iceland  and  Azores  plumes  in  the  North  Atlantic,  and  the  Tristan  plume  in  the  south 
Atlantic.  Such  a  scenario  implies  that  plumes  are  a  major  source  of  lithospheric  accretion  as 
proposed  by  Morgan  and  Smith  [1992]  and  Morgan  et  al.  [1995], 

CONCLUSIONS 

In  our  numerical  investigations  of  steady-state  stationary  ridges,  we  have  derived 
scaling  laws  consistent  with  those  of  [Ribe,  1996],  indicating  that  they  are  insensitive  to 
differences  in  numerical  method  or  model  boundary  conditions.  Plume  width  W  and 
maximum  plume  ridge  communication  distance  increase  with  the  plume  width  scale 
and  modified  plume  buoyancy  number  Tlby.  In  the  case  of  a  migrating  ridge,  the 
distance  of  plume-ridge  interaction  is  reduced  when  the  ridge  migrates  toward  the  plume 
due  to  the  excess  drag  of  the  leading  plate.  After  the  ridge  passes  over  and  migrates  away 
from  plume,  the  distance  of  plume-ridge  interaction  is  enhanced  due  primarily  to  the 
reduced  drag  of  the  slower-moving  trailing  plate,  and  secondarily  to  the  pattern  of  thermal 
erosion  of  the  lithosphere. 

To  test  our  plume-ridge  models  we  compare  model  predictions  of  along-isochron 
mantle-Bouguer  and  bathymetric  anomalies  with  observations  of  the  Galapagos  plume- 
migrating  ridge  system.  The  models  predict  the  amplitudes  and  widths  of  the  observed 
anomalies  with  a  plume  source  temperature  anomaly  of  80-120°C,  radius  of  80-100  km, 
and  volume  flux  of  4.5x10^  km^/m.y.  The  models  also  predict  the  approximate  increase  in 
anomaly  amplitudes  with  isochron  age  which  reflects  increased  crustal  production  in  the 
past  when  the  ridge  was  closer  to  the  Galapagos  plume.  Crustal  production  rates  of  the 
Galapagos  Islands,  as  estimated  from  the  observed  island  topography,  are  also  matched 
reasonably  well  by  model  predictions.  Predictions  of  the  geochemical  signature  of  the 
plume  along  the  present-day  ridge  suggest  that  mixing  between  the  plume  OIB  and  ambient 
MORE  source  does  not  occur  in  the  asthenosphere  but  instead  most  likely  occurs  deeper, 
possibly  by  entrainment  of  the  depleted  mantle  as  the  plume  ascends  from  its  deep  source 
region.  These  numerical  models  suggest  that  plumes  may  spread  perpendicular  to  the 
direction  of  plate  motion  over  distances  125-200%  broader  than  the  maximum  distance  to 
which  they  interact  with  ridges.  Plumes  may  therefore  comprise  a  significant  percentage  of 
the  oceanic  lithosphere. 
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Table  1.  Notation 


Variable 

Meaning 

Value 

Units 

specific  heat 

1000 

J  kg-'  °c-' 

D 

fluid  depth 

400 

km 

E 

activation  energy 

1.9x105 

J 

g 

acceleration  of  gravity 

9.8 

m/s^ 

M 

melt  fraction 

wt% 

P 

pressure 

Pa 

P 

plume  tracer  concentration 

Q 

volumetric  plume  flux 

km^/m.y. 

Qr 

fraction  of  plume  flux  crossing  the  ridge 

R 

gas  constant 

8.314 

J  K-'  mol-' 

So 

characteristic  plume  thickness 

(48(2Po/gAp)>/4 

km 

AS 

entropy  change  on  melting 

400 

J  kg-‘  °C 

T 

mantle  potential  temperature 

°C 

Tr 

mantle  real  temperature 

K 

plume  temperature  anomaly 

°c 

u  (u,v,w) 

mantle  flow  rate  vector 

km/m.y. 

U 

ridge  full  spreading  rate 

km/m.y. 

V 

activation  volume 

4x10-6 

m3 

w 

along-axis  plume  width 

km 

Wo 

characteristic  plume  width 

{QlW^ 

km 

Xp 

plume-ridge  distance 

km 

^max 

maximum  distance  of  plume-ridge  interaction 

km 

X 

melt  depletion 

wt% 

a 

coefficient  of  thermal  expansion 

3.4x10-5 

K-^ 

p 

coefficient  of  depletion  density  reduction 

0.024 

Y 

r]olnp 

K 

thermal  diffusivity 

31 

kmVm.y. 

n 

viscosity 

Pa  s 

Vo 

reference  viscosity 

Pa  s 

Vp 

plume  viscosity  at  0.5Z) 

Pa  s 

Hh 

buoyancy  number 

Qallfl 

Ho 

upslope  number 

gm^yVi^miu 

P 

mantle  density 

kg/m^ 

Pmax 

ridge  crustal  density  closest  to  the  plume 

2900,  3000 

kg/m^ 

Pm 

melt  density 

2900 

kg/m^ 

Po 

mantle  reference  density 

3300 

kg/m^ 

a 

buoyancy  scaling  parameter 

g  Apt  Air]  0 

1/ms 
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Table  2.  Experimental  parameters  and  scaling  quantities 


Run 

W  /  U  (km/my) 

ATpCC) 

r 

Q’/Q{\0^  km-Vmy) 

W'Up=0)  /  W(;c^=0)  (km) 

1 

516/40 

100 

1.0 

0.61 

55  /  0.68 

0.33  /  131 

0.63  /  250 

2 

774  /  60 

100 

1.0 

0.29 

58  /  0.72 

0.27  /  110 

0.56  /  225 

3 

1290  /  100 

100 

1.0 

0.12 

64  /  0.79 

0.22/  89 

0.50  /  200 

4 

386  /  30 

200 

1.0 

3.97 

100  /  1.24 

0.51  /  203 

1.00  /  400 

5 

744/  60 

200 

1.0 

1.04 

105  /  1.30 

0.37  /  147 

0.75  /  300 

6 

1290  /  100 

200 

1.0 

0.40 

112  /  1.39 

0.29  /  118 

0,63  /  250 

7 

516  /  40 

100 

2.352 

0.60 

54  /  0.67 

0.32  /  129 

0.69  /  275 

8 

774  /  60 

100 

2.352 

0.28 

57  /  0.70 

0.27  /  108 

0.63  /  250 

9 

1290  /  100 

100 

2.352 

0.1 1 

61  /  0.75 

0.22  /  87 

0.50  /  200 

10 

386  /  30 

200 

5.053 

3.13 

157  /  1.95 

0,64  /  255 

1.65  /  660 

1  1 

774  /  60 

200 

5.053 

1.25 

126  /  1.56 

0.40  /  161 

0.94  /  375 

12 

1290  /  100 

200 

5.053 

0.46 

130  /  1.61 

0.32  /  127 

0.75  /  300 

Primes  denote  dimensionless  quantities  and  are  listed  adjacent  to  their  scaled  quantities.  Input  parameters  are  f/, 
ATp,  and  rheology  law,  which  controlled  7.  The  remaining  quantities  are  model  output.  Runs  4  and  10  had 
numerical  box  dimensions  3.2  D  x  2.0  D  x  1.0  £)  with  128  x  64  x  50  grids  in  x,  y,  and  z  respectively.  The  other 
runs  had  box  dimensions  3.2  D  x  1.0  £>  x  \.0  D  with  128  x  32  x  50  grids  in  x,  y,  and  z  respectively.  Rayleigh 
number  was  1.83  x  10^  based  on  T„=  1300  °C  and  r}„  =  5  x  10^^  Pa  s. 
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Figure  1 .  Perspective  diagram  illustrating  steady-state  flow  (small  arrows)  and  potential 
temperature  fields  (contoured  at  50°C-intervals  for  r>  1300°C  and  100°C  intervals  for  r< 
1300  °C)  of  an  example  calculation  with  ATp  =  100°C  and  t/  =  60  km/m.y.  (experiment  7, 
Table  2).  The  ridge  axis  is  located  at  x  =  320  km,  the  plume  source  is  centered  at  x  =  450 
km.  The  maximum  vertical  velocity  is  1 15  km/my.  Both  top  (z  =  0)  and  bottom  (z  =  D) 
boundaries  are  isothermal  planes  with  the  bottom,  a  free  slip  boundary  and  the  top,  fixed  at 
a  horizontal  velocity  of  UH  (large  horizontal  arrow)  at  x  >  320  km  and  -f//2  at  x  <  320  km. 
All  boundaries  are  closed  to  flow  both  in  and  out  of  the  numerical  box,  thus  material  flows 
downward  at  the  ends  of  the  box  and  recirculates  toward  the  ridge  axis  along  the  base  of 
the  box.  The  effect  of  this  recirculation  on  the  interaction  between  plume  and  ridge  are 
insignificant.  Note  the  cooling  lithosphere  which  slopes  towards  the  ridge  axis. 
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Figure  2.  Model  predictions  of  scaled  along-axis  width  versus  modified  buoyancy  number 
TJhY.  Open  circles  are  for  runs  with  7=  1.0  and  temperature  anomalies  100  and  200°C. 
Gray  circles  are  for  fully  pressure-  and  temperature-dependent  plume  viscosity  calculations 
with  7  =  2.35  and  ATp  =  100°C,  and  black  circles  are  for  temperature-dependent  plume 
viscosity  calculations  with  7  =  5.05  and  ATp  =  200°C.  The  curve  is  the  best  fitting  scaling 
law  described  by  Eq.  6. 
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Figure  2 
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Figure  3.  Steady  state  isosurface  of  plume  tracer,  P  =  0.4,  as  viewed  from  the  bottom  of 
the  box  looking  upward  (shading  denotes  illumination  from  the  right  of  the  figure).  Small 
arrows  illustrate  horizontal  velocities  in  the  horizontal  plane  at  z  =  64  km.  The  ridge  axis  is 
marked  by  the  bold  line.  Experimental  conditions  are  those  of  experiment  7,  the  same  as  in 
Fig.  1.  a)  =  0,  b)  Xp  =  100,  c)  Xp  =  150.  Note  that  the  width  of  the  plume  at  the  ridge 
axis  decreases  with  increasing  Xp.  The  maximum  distance  to  which  the  plume  reaches  the 
ridge  is  Xp  =  150. 


AIong-Axis  (km) 


Figure  4.  Numerical  results  of  along-axis  plume  width  (scaled  by  width  for  Xp  =  0)  versus 
scaled  plume-ridge  distance.  As  in  Fig.  2  open  circles  are  for  runs  with  7=  1.0,  gray 
circles  are  for  7=  2.35,  black  circles  are  for  7=  5.05.  (a)  The  best  fitting  polynomials  of 
the  form  given  in  Eq.  7  are  shown  for  7=  1.0  (solid),  7=  2.35  (dashed),  and  7=  5.05 
(dotted).  The  different  widths  of  the  curves  illustrate  the  dependence  of  F2  on  7  (b)  Same 
as  in  (a)  but  including  F2  (Eq.  8)  which  collapses  the  points  onto  to  a  single  curve.  The 
standard  deviation  misfit  of  Eq  8  to  the  numerical  results  is  0.13.  The  mismatch  to  the 
numerical  points  for  Xp/(WoF2)  <  ~1.0  may  suggest  a  dependence  on  higher  order  terms  of 
Xp  which  we  chose  not  attempt  to  resolve. 
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W/W(Xp=0)  S  W/W(Xp==0) 


Figure  4 


Figure  5.  Dependence  of  the  ratio  of  plume  volume  flux  crossing  the  ridge  Qr  versus 
scaled  plume  ridge  distance  divided  by  stretching  function  ¥2-  Symbols  are  the  same  as  in 
Fig.  2.  The  solid  line  is  the  best  fitting  line  of  Eq.  9  with  a  standard  deviation  misfit  of 
0.08.  The  mismatch  to  the  numerical  points  for  XpHy^ ^  ~0.7  may  suggest  a 
dependence  on  higher  order  terms  of  Xp  which  we  chose  not  attempt  to  resolve. 
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Figure  6.  Contours  of  lithospheric  erosional  thickness  are  normalized  by  characteristic 
plume  thickness  =  128  km  for  experiment  7  (same  as  Fig  3c).  The  ridge  axis  is 

marked  by  the  shaded  vertical  line  and  the  plume  source  is  shown  as  the  gray  semicircle  at 
Xp  =  1 50.  The  thickness  of  the  lithospheric  rheological  boundary  layer  is  defined  as  the 
depths  over  which  77/77  0  ^  10-  Erosional  thickness  is  the  difference  between  the  boundary- 
layer  thickness  above  the  plume  and  that  of  normal  lithosphere  as  defined  along  the  ridge- 
perpendicular  profile  at  y  =  l.OZ). 


130 


Across-Axis  (km) 


-100  0  100  200  300  400 


-10  12  3  4 


Figure  6 


131 


Along-Axis  (km) 


Figure  7.  Temperature  fields  (contoured  at  50  °C  intervals  in  the  plume  and  at  100°C 
intervals  in  the  lithosphere)  and  velocities  (arrows)  in  across-axis,  depth  cross-sections 
through  the  center  of  the  plume  source  (y  =  0).  Experimental  parameters  are  the  same  as  in 
Fig  1 .  (experiment  7)  but  the  ridge  is  migrating  in  the  positive  x^-direction  at  the  half 
spreading  rate  of  30  km/m.y.  (a)  Ridge  is  migrating  toward  the  plume  therefore  the  plume 
is  beneath  the  faster  moving  Plate  1.  (b)  Ridge  is  migrating  away  from  the  plume  therefore 
the  plume  is  beneath  the  stationary  Plate  2. 
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Figure  8.  Numerical  results  of  scaled  along- axis  width  versus  scaled  plume-ridge  distance 
for  migrating  ridge  cases  of  experiment  7  with  f//2=  30  km/m.y.  The  bold  curve  as 
defined  in  Eq.  7  is  that  predicted  for  steady  state  conditions  with  a  stationary  ridge.  Open 
triangles  are  for  =  10  km/m.y.  shown  with  best  fitting  (solid)  curve  of  the  form  in  Eq. 
11;  gray  triangles  are  for  Vr  =  20  km/m.y.  shown  with  best  fitting  (dotted)  curve;  solid 
triangles  are  for  Vy  -  30  km/m.y.  shown  with  best  fitting  (dashed)  curve.  Mismatches  are 
largest  near  the  apexes  of  the  curves  and  are  due  in  part  to  difficulty  in  resolving  curvature 
where  slope  in  W  is  small,  and  to  a  possible  dependence  on  higher  order  terms  of  Xp  which 
we  chose  not  consider. 
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Figure  8 


Figure  9.  (a)  Numerical  results  showing  the  dependence  of  F5  on  Uhyand  scaled  ridge 
migration  rate.  Open  circles  are  for  runs  with  7=  1.0,  gray  circles  are  for  7=  2.35,  black 
circles  are  for  7=  5.05.  The  line  is  the  best  fitting  function  of  Eq.  12.  (b)  Numerical 
results  showing  the  dependence  of  Fg  on  Ubyand  scaled  ridge  migration  rate.  Circles  are 
patterned  as  in  (a).  The  curve  is  the  best  fitting  function  of  Eq.  13. 
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Figure  10.  Contours  of  lithospheric  erosion  thickness  normalized  by  characteristic  plume 
thickness  =  120  km  for  the  migrating  ridge  cases  of  experiment  7  with  U  =  30 

km/m.y.  The  plume  source  (shaded)  is  now  at  jc^  =  -170  km  beneath  the  slower  moving 
Plate  2.  (a)  The  region  of  erosion  is  broadest  for  the  case  where  y^=  10  km/m.y.  The  area 
of  erosion  becomes  more  confined  to  the  plume  source  with  increasing  ridge  migration 
rates  (b)  VV=  20  km/m.y.  and  (c)  Vr=  30  km/m.y. 
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Along-Axis  (km) 


Figure  1 1 .  Map  of  the  regional  bathymetry  of  the  Galapagos  plume-ridge  system 
(shipboard  and  EtopoS  bathymetry  from  Ito  and  Lin  [1995a]).  The  present-day  ridge  axis 
is  the  southern-most  set  of  white  lines,  and  the  isochrons  of  Ito  and  Lin  [1995a]  (taken 
from  Wilson  and  Hey  [1995])  are  shown  to  the  north  on  the  Cocos  Plate.  The  plume 
center  is  taken  to  be  the  eastern-most  island  Fernandina  as  shown  by  the  circle  of  radius 
100  km.  The  dashed  box  shows  the  region  of  bathymetry  used  to  calculate  the  crustal 
volume  flux  of  the  archipelago. 
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Figure  12.  Comparisons  between  observed  (thick  gray)  and  predicted  along-isochron, 
bathymetric  profiles  from  model  1  (solid)  and  model  2  (dashed).  Model  profiles  are  the 
combined  isostatic  topography  of  axial  crustal  thickness  and  mantle  density  variations. 
Maximum  values  of  crustal  thickness  predicted  by  models  1  and  2  are  labeled  as  ACri  2- 
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Figure  13.  Comparisons  between  observed  (thick  gray)  and  predicted  along-isochron, 
MBA  profiles  from  model  1  (solid)  and  model  2  (dashed). 
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Figure  14.  Total  amplitude  of  along-isochron  (a)  bathymetry  and  (b)  MBA  variations  are 
plotted  versus  isochron  age.  Thick  gray  lines  are  observed  variations,  solid  lines  are 
variations  predicted  by  model  1,  and  dashed  lines  are  variations  predicted  by  model  2.  The 
pairs  of  model  curves  are  those  assuming  a  crustal  density  of  2900  kg/m^  (upper-bound) 
and  3000  kg/m^  (lower-bound)  at  the  point  of  the  ridge  closest  to  the  plume  (91°W).  (c) 
Crustal  volume  flux  of  the  Galapagos  Archipelago  versus  age  as  predicted  from  model  1 
(solid)  and  model  2  (dashed)  is  compared  with  crustal  volume  fluxes  calculated  by 
assuming  isostatic  compensation  of  the  island  topography  (solid  gray)  (see  text).  The 
upper-bound  gray  curve  is  that  assuming  a  crustal  density  of  3000  kg/m^  and  the  lower- 
bound  curve  is  that  assuming  a  crustal  density  of  2900  kg/m^. 
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Figure  15.  Observed  variations  (dots)  in  [La/Sm]ef  from  ^c/z/Z/mg  et  al.  [1982]  is 
compared  with  accumulated  concentration  of  plume  tracer  along  the  ridge  (Eq.  19)  for 
model  1  (solid)  and  model  2  (gray). 
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Conclusions 


Bathymetric  and  gravity  observations  at  five  prominent  plume-ridge  systems  reveal 
broad  wavelength  anomalies  that  reflect  anomalously  low  density  subsurface  structure 
imposed  by  the  near-ridge  plumes.  Along-axis  bathymetry  shallows  by  as  much  as  4.5  km 
toward  the  plumes  while  along-axis  mantle-Bouguer  anomalies  become  increasingly 
negative  by  as  much  as  -300  mGal  toward  the  plumes.  We  estimate  that  -70%  of  the 
anomaly  amplitudes  are  due  to  thickened  axial  crust  with  the  remaining  -30%  due  to 
anomalously  low  density  mantle,  both  of  which  are  caused  by  anomalously  hot  mantle 
temperatures  imposed  by  the  near-ridge  plumes. 

The  amplitudes  of  along-isochron  MBA  and  bathymetric  anomalies  are  largest  at  the 
ridge-centered  Iceland-Mid- Atlantic  Ridge  (MAR)  system  and  at  the  ridge-centered  cases  of 
the  Tristan-MAR  system.  The  anomaly  amplitudes  decrease  with  increasing  plume-ridge 
distance  most  likely  reflecting  reduced  crustal  production  as  the  ridges  migrated  away  from 
the  plumes.  At  a  plume-ridge  distance  of  -500  km  the  available  data  at  the  Tristan-MAR 
system  show  no  discernable  anomalies  suggesting  a  maximum  distance  that  these  plumes 
affect  ridge  structure  significantly.  Residual  bathymetric  anomaly  widths  along  the 
isochrons,  however,  appear  to  be  most  sensitive  to  spreading  rate  and  decrease  with 
increasing  spreading  rate  from  2700  km  at  the  slow  spreading  Iceland-MAR  system,  to  < 
500  km  at  the  fast  spreading  Easter-EPR  system. 

While  the  above  studies  place  constraints  on  the  amplitude  and  extent  of  plume-imposed 
subsurface  density  anomalies,  our  numerical  modeling  studies  examine  the  possible  causes 
of  such  anomalies.  Numerical  models  of  ridge-centered  plumes  indicate  that  along-axis 
plume  width  W  scales  with  plume  volume  flux  Q,  ridge  full  spreading  rate  U, 
ambient/plume  viscosity  ratio  7,  and  buoyancy  number  Ut  according  to  W  = 
1.2,l{QIU)^'Hnb  7,)®-®'^.  Thermal  buoyancy  is  the  most  important  driving  force  while 
melting  effects  of  latent  heat  loss,  depletion  buoyancy,  and  melt-retention  buoyancy  yield 
competing  effects  which  do  not  change  the  above  scaling  argument.  Numerical  simulations 
of  the  Iceland-MAR  system  suggest  two  end-member  source  models.  The  first  model  has  a 
source  radius  of  300  km,  temperature  anomaly  of  75°C,  and  volume  flux  of  1.2  x  10^ 
km^/m.y.,  while  the  second  has  a  source  radius  of  60  km,  temperature  anomaly  of  170°C, 
and  volume  flux  of  2.1  x  10^  km^/m.y.  The  first  model  explains  well  the  observed  crustal 
thickness,  bathymetric,  and  MBA  variations  along  the  MAR  and  Iceland,  but  the  second 
model  requires  substantial  along-axis  melt  transport  in  order  to  explain  the  observations. 
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This  second  model  may  be  more  representative  of  the  Iceland  plume  based  on  similarities 
between  predicted  and  observed  mantle  P-wave  anomalies. 

For  off-axis  plumes,  along-axis  width  scales  again  with  and  Ut,  7,  in  a  similar 

manner  to  the  ridge-centered  plume  case.  For  steady-state  plumes  near  stationary  midocean 
ridges  W  decreases  with  increasing  plume-ridge  distance  and  becomes  zero  at  a  maximum 
plume-ridge  interaction  distance  x^ax^  which  increases  with  and  17^  7.  When 

ridges  migrate  toward  plumes,  however,  predicted  values  of  W  and  x^ax  are  reduced 
relative  to  the  case  of  a  stationary  ridge  by  as  -24%  due  to  the  enhanced  drag  of  the 
overlying  plate  that  leads  the  migrating  ridge.  On  the  other  hand,  when  ridges  migrate 
away  from  plumes,  W  and  x^ax  are  predicted  to  increase  by  -36%  due  to  reduced  drag  of 
the  trailing  plate;  enhanced  erosion  of  the  lithosphere  also  enhances  W  and  Xmax  but  to  a 
degree  that  is  secondary  to  the  effects  of  reduced  plate  shearing.  Numerical  models  of  the 
Galapagos  plume-ridge  system  predict  MBA  and  bathymetric  anomalies  that  match 
successfully  the  amplitudes  and  widths  of  the  observed  anomalies,  as  well  as  the  increase 
in  anomaly  amplitude  with  isochron  age.  The  implied  Galapagos  plume  source  has  of 
radius  80-100  km  and  temperature  anomaly  of  80-1 20°C.  In  addition,  predicted  chemical 
signatures  of  the  plume  along  the  model  ridge  suggest  that  mixing  between  the  plume  and 
ambient  mantle  occurs  deeper  than  the  asthenosphere,  most  likely  due  to  entrainment  of  the 
ambient  mantle  as  the  plume  ascends  from  its  deep  source  reservoir. 

Thus  for  a  few  prominent  plume-ridge  systems,  we  have  begun  to  quantify  the 
influence  of  near-ridge  plumes  on  ridge  crustal  and  mantle  density  structure.  The 
suggestion  that  subsurface  structure  along  seafloor  isochrons  reflects  past  interaction 
between  plumes  and  ridges  warrants  further  investigations  to  test,  in  the  form  of  land-based 
data  analyses  as  well  as  sea  going  geophysical  and  geochemical  surveys.  In  addition,  we 
have  learned  a  great  deal  about  how  mantle  flow  might  behave  at  plume-ridge  systems. 
Our  models  also  require  further  studies  to  test,  most  likely  with  mantle  seismological 
studies.  Finally,  as  we  have  discussed  for  the  Iceland-MAR  system,  along-axis  melt 
transport  may  be  a  first-order  process  for  this  and  possibly  other  plume-ridge  systems.  If 
this  is  the  case,  it  may  be  time  to  re-examine  our  established  ideas  of  plume-ridge 
interaction  and  possibly  crustal  accretionary  processes  at  ridges  in  general. 


Appendix 

Laboratory  Investigation  of  the  Interaction  of  Off-Axis  Mantle 

Plumes  and  Spreading  Centres 
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Mantle  plumes  and  mid-ocean  ridge  spreading  centres  are  the 
dominant  phenomena  through  which  mass  and  heat  are  transported 
from  the  mantle  to  the  Earth’s  surface.  It  now  seems  that  the 
dispersion  of  near-ridge  plumes  beneath  the  lithosphere  is  modu¬ 
lated  strongly  by  mid-ocean  ridges;  in  particular,  geochemical  and 
geophysical  ob^rvations  have  suggested  that  rising  plumes  are 
diverted  towards  and  feed  nearby  ridges*'^.  Here  we  confirm  the 
feasibility  of  this  model  with  laboratory  experiments  that  incorpor¬ 


ate  the  essential  physical  and  fluid  dynamic  aspects  of  a  plume- 
ridge  upper  mantle  system.  Our  results  indicate  that  an  off-axis 
plume  may  communicate  thermally  and  chemically  with  a  spread¬ 
ing  ridge  through  a  narrow,  sub-horizontal  conduit  instead  of  a 
broader,  radially  spreading  plume  head.  A  necessary  condition  for 
this  communication  is  the  presence  of  a  lithospheric  or  rheological 
boundary  layer  that  thickens  away  from  the  ridge  axis  owing  to 
conductive  cooling.  Interestingly,  we  find  that  for  high  plume  tem¬ 
peratures,  increasing  the  plume  thermal  buoyancy  may  inhibit 
rather  than  enhance  plume-ridge  interaction,  as  a  result  of 
increased  erosion  of  the  overlying  lithosphere. 

Recent  laboratory®  and  numerical  experiments^  have  consid¬ 
ered  the  dynamics  of  plume-ridge  interaction  for  the  ridge- 
centred  case;  however,  the  difficult  question  remains  as  to  how 
a  plume  and  a  ridge  interact  thermally,  chemically  and  dynami¬ 
cally  when  the  plume  is  located  off  axis.  A  model  of  sub-horizon¬ 
tal  pipe-like  flow  from  the  off-axis  plumes  to  a  ridge^  axis  along 
the  base  of  the  rigid  lithosphere  has  been  suggested*  (Fig.  \a). 
Geochemical  studies*®'"^  and  two-dimensional  numerical  experi¬ 
ments’’"*^  support  this  channel-flow  model  and  suggest  that  geo¬ 
chemical  communication  may  persist  over  long  periods  of  time 
and  plume-ridge  separation  distances  as  high  as  1,200  km  (ref. 
16).  This  laboratory  experimental  project  is  the  first  fully  three- 


D.C.  drive 


FIG.  1  a.  This  diagram  illustrates  the  conceptual  model  that  near¬ 
ridge  plumes  rise  and  then  flow  toward  ridges  along  the  base  of  the 
RBL  from  A  to  S  (ref.  16).  Up  is  vertical  velocity  within  the  plume 
conduit,  b.  Diagram  of  experimental  apparatus.  Mylar  sheeting  is 
pulled  along  the  fluid  surface  to  simulate  plate  spreading.  From 
source  reels,  the  Mylar  is  threaded  through  two  bars  at  the  spreading 
axis  and  at  the  tank  edges  (take-up  bars).to  ensure  contact  between 
sheeting  and  working  fluid.  Mylar  is  then  pulled  around  take-up  bars 
by  a  synchronous  high-torque  d.c,  motor.  The  fluid  surface  is  cooled 
by  circulating  fluid  from  a  refrigerated  cold  bath  through  a  series  of 
(70  cm  X 10  cm  X  2  cm)  metal  jackets  suspended  2  mm  above  the 
fluid.  Tank  sidewalls  and  fluid  surface  are  insulated.  Plume  flow  is 
monitored  with  shadowgraph  and  time-lapse  laser  photographs  from 
ridge-perpendicular  and  ridge- pa  rail  el  viewpoints.  Working  fluid 
viscosity  follows  an  Arrhenius-type  law  of  the  form  fi  ~  exp  ((1,888/ 
(T+  93.3))  - 11.48],  where  T  is  temperature  in  ”0  and  /i  is  viscosity 
(in  Pa  s).  c,  This  close-up  slice  through  the  tank  centre  illustrates  the 
configuration  of  the  ridge  axis  and  plume  source.  Arrows  illustrate 
hypothetical  fluid  flow.  Nine  RTDs  are  positioned  at  5-cm  Intervals 
along  the  ridge  axis  to  monitor  axial  temperature  variability  (note 
that  ridge  RTD  5  lies  on  an  orthogonal  line  from  the  plume  source). 
Before  running  experiments,  we  allow  the  fluid  to  equilibrate  for 
several  days  at  room  temperature,  —20  °C.  We  then  establish  large- 
scale  plate-driven  flow  by  running  the  Mylar  drives  for  60  min  before 
heating  the  plume  source.  In  experiments  2  and  4,  surface  coolers 
are  maintained  at  10  °C  for  60  min  to  produce  an  RBL  before  activ¬ 
ating  the  Mylan  in  experiment  3,  coolers  are  maintained  at  0  *0 
for  120  min  before  activating  the  Mylar.  Experiments  run  for  150- 
280  min,  depending  on  whether  or  not  there  is  surface  cooling. 
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FIG.  2  Shadowgraph  photographs  at  three  instants 
during  experiment  1.  Variations  in  fluid-temperature 
gradients  focus  the  illuminating  light  to  yield  a  bright 
halo  at  the  top  of  the  plume  head  (beneath  black 
curves),  a.  31  min  after  turning  on  the  disc  heater, 
the  plume  has  separated  from  the  source  at  the 
base  of  the  tank,  generating  a  single  plume  charac¬ 
terized  by  a  broad  leading  head  (-6  cm  across)  and 
a  narrow  trailing  conduit  (-1.8  cm  wide),  b,  At 
41  min.  the  plume  begins  impinging  upon  the  fluid 
surface.  The  maximum  ridgeward  deflection  due  to 
the  plate-driven  return  flow  is  -2.2  cm,  or  roughly  a 
conduit  width.  The  ascent  rate  of  the  plume  head  is 
-0.8  cm  min"\  or  close  to  twice  the  plate  speed. 
Fluid  velocity  in  the  conduit  exceeds  3  cm  min'^  as 
measured  by  tracking  neutrally  buoyant  Delrin  beads 
(not  shown)  released  periodically  at  the  plume 
source,  c,  At  72  min.  the  plume  head  has  stalled 
and  flattened:  head  and  conduit  are  being  sheared 
away  from  the  ridge  near  the  fluid  surface. 


Plume  Heat  Source 


dimensional  variable-viscosity  study  of  this  problem,  and  it 
exposes  the  mechanisms  by  which  an  off-axis  plume  overcomes 
the  lithospheric  drag  that  draws  material  away  from  a  ridge,  to 
successfully  feed  the  nearby  spreading  centre. 

To  test  the  conceptual  model  of  plume-ridge  channelling  (Fig. 
la),  we  simulate  a  plume-ridge,  upper  mantle  system  with  a  tank 
of  a  concentrated  sucrose  solution  which,  like  the  Earth’s  mantle, 
is  strongly  temperature  dependent  (Fig.  \b).  Plate-driven  mantle 
Row  is  simulated  by  dragging  two  Mylar  sheets  in  opposite 
directions  on  the  surface  of  the  fluid  at  a  constant  rate  of 
0.35  cm  min“’  (Uo,  half-spreading  rate).  Buoyantly  driven  flow 
is  produced  through  a  supply  of  thermal  energy  from  a  disc¬ 
shaped  resistance  heater  (that  is,  the  plume  source)  positioned 
at  the  base  of  the  tank.  The  two  parameters  we  vary  are  the 


surface  temperature,  thus  the  thickness/age  of  the  upper  rheo¬ 
logical  boundary  layer  (RBL),  and  the  plume  source  tempera¬ 
ture,  controlling  the  strength  of  the  rising  plume.  Fluid 
temperature  is  continuously  monitored  at  the  disc  heater,  at  the 
surface  of  the  Mylar,  along  a  vertical  profile  within  the  fluid, 
and  along  the  ridge  axis  with  resistance  temperature  detectors 
(RTDs)  (Fig.  Ic). 

Plume  buoyancy  is  caused  by  density  reduction  of  the  plume 
fluid  due  to  thermal  expansion  according  to 

Ap=={Po~pp)  =  apo{Tp-To)  (1) 

where  or  is  the  expansion  coefficient  (4.6  x  10“"^),  and  Po,  To 
and  pp ,  Tp  are  reference  and  plume  densities  and  temperatures. 


TABLE  1  Plume-ridge  experiments 


(a)  Parameters 


Surface  temp. 

Plume  source 

Plume  density 

Mean  plume 

Plume 

Plume  buoy. 

contrast  from 

temp,  contrast 

contrast. 

conduit 

buoyancy 

flux  to  ridge*, 

ambient  of  20  X 

RBL  slope 

from  ambient 

Ap/po 

viscosity,  pp 

number. 

B./8s 

Exp.  no. 

(^'C) 

(deg) 

(X) 

(%) 

(Pas) 

Bn 

(%) 

1 

0 

0 

40 

1.8 

4.9 

52 

0 

2 

-6 

1 

40 

1.8 

4.9 

52 

3 

3 

-10 

3 

40 

1.8 

4.9 

52 

10 

4 

-6 

1 

48 

2.2 

2.5 

60 

0 

(5)  Comparison  of  laboratory  and  expected  mantle  parameter  ranges 

Pu/Po 

Po/Pd  Ap/po  (%) 

Bn 

Pep 

Pep/Pe 

Laboratory 

10 

50-100 

1.8-2.2 

50-60 

1.750-2,500 

10 

Mantle 

>10^ 

-100  -0.5-1.5 

5-60 

10^-10^ 

1-100 

*  Calculated  as  a  ratio  of  plume  buoyancy  flux  at  the  ridge,  B,-  poaUoDrb\  f  (T(x)-To)  dx,  where  the  x  axis  is  along  the  ridge,  0^,  is  mean  RBL 
thickness  and  T{x)  is  along-axis  temperature,  and  plume  source  buoyancy  flux^^,  8^  =  poaATlJp7rid/2f,  where  Up  is  bead  velocity  within  the  plume 
conduit  and  d  is  conduit  diameter  (1.5-1.8  cm).  8s  ranges  between  0.003-0.004  gs  ^.  Reference  density,  po,  at  ambient  temperature  (20  X)  is 
1.4  g  cm  ^  Important  parameters  for  comparing  results  on  8, /8s  are  highlighted  in  bold.  Comparisons  are  also  made  using  plume  Peclet  number, 
Pep-UpD/K,  and  the  ratio  Pep/Pe.  Mantle  PCp  is  calculated  from  Whitehead  and  Luther's^^  equation  for  Up.  Bn,  Pep  and  PCp/Pe,  which  best 
represent  tl^  vigour  of  plume  convection  relative  to  plate-driven  flow,  and  Po/Pp  scale  well  with  expected  mantle  values,  pl/p©  and  Ap/p©  reflect 
differences4n  laboratory  and  mantle  material  properties. 
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FIG.  3  a.  Photograph  of  the  tank  fluid  during  experiment  2  showing 
locations  of  7  Defrin  beads  (trajectories  sketched  in  black).  Bead  no.  1 
reached  the  high -viscosity  upper  boundary  layer  then  travelled  along  a 
sub-horizontal  path  to  the  ridge  axis  similar  to  the  depiction  in  Fig.  la. 
This  photograph  was  taken  65  min  after  plume  initiation,  50  min  after 
the  plume-RBL  impact,  and  15  min  after  bead  no.  1  hit  the  ridge.  The 
bead  is  now  frozen  into  the  migrating  plate.  The  photograph  is  taken 
from  a  mid-depth  fluid  level,  and  the  dark  sloping  line  is  the  ridge  axis 
(line  of  RTDs)  as  viewed  from  below,  through  the  fluid.  The  bead  below 
and  to  the  left  of  no.  7  was  introduced  while  setting  the  bead  source 
and  is  not  part  of  the  experiment  b.  Fluid  temperatures  measured  by 
eight  RTDs  along  the  ridge  axis  at  different  times  during  experiment  2, 
showing  evolution  of  a  narrow  axial  anomaly.  Temperatures  increase 
with  time  primarily  at  RTD  5,  peaking  at  -3  “C  above  ambient  tempera¬ 
ture.  Temperatures  are  still  increasing  at  95  min,  or  45  min  beyond 
bead  no.  I’s  arrival,  c,  Fluid  temperatures  along  the  ridge  axis  at  the 
conclusions  of  the  four  experiments.  Without  the  RBL  (experiment  1)  or 
if  the  plume  is  too  hot  (experiment  4),  the  plume  fails  to  reach  the  ridge 
and  enhance  ridge  temperatures.  Note  the  broader  axial  anomaly  for 
the  case  with  a  larger  (3°)  RBL  slope  (experiment  3). 


respectively.  We  characterize  the  strength  of  the  plume  by  the 
dimensionless  buoyancy  number®,  ,  defined  as 


(2) 


where  ^  is  the  acceleration  due  to  gravity,  Po  is  reference  dynamic 
viscosity  and  Q  is  the  volumetric  plume  fiux“°,  which  is  a  con¬ 
stant  0.14  cm^s"'  between  experiments.  Characterizing  plumes 
in  this  manner  enables  us  to  compare  plumes  quantitatively 
between  experiments  (see  Table  1)  and  provides  a  measure  of 
how  well  our  laboratory  plumes  represent  Earth  examples.  Our 
laboratory  Bn  are  near  the  upper  limit  of  the  expected  range  for 
the  Earth  of  7-59  (ref.  8). 


Length  and  time  scales  in  the  laboratory  models  are  related 
to  the  mantle  through  the  Peclet  number, 

?c=UoD/k  (3) 

Thermal  diffusivity,  k,  of  the  laboratory  fluid  is  0.001  cm^  s“’ 
and  the  corresponding  mantle  value  is  0.01  ernes'"'.  We  define 
the  length  scale,  D,  as  the  thickness  of  the  laboratory  fluid 
( 1 7  cm)  corresponding  to  600  km  of  the  upper  mantle.  Thus, 
our  0.35  cm  min”’  Mylar  speed  yields  Pe=  100  and  scales  to  a 
slow  mantle  full-spreading  rate  of  ->>'1  cm  yr”'.  Likewise,  our 
laboratory  plume-ridge  separation  distance  of  14  cm  scales  to  a 
mantle  distance  of  '^500  km. 

The  four  experiments  we  present  here  are  selected  from  a  total 
of  nine  to  highlight  the  relative  roles  of  surface  cooling  (compare 
experiment  1  with  2,  and  2  with  3)  and  plume  source  temperature 
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(compare  experiments  2  and  4)  on  plume-ridge  communication 
(see  Table  1).  In  experiment  1  (Fig.  2),  the  fluid  surface  is  not 
cooled,  but  rather  maintained  at  room  temperature;  therefore 
no  sloping  RBL  is  present.  As  the  rising  plume  reaches  the  sur¬ 
face,  the  plume  head  flattens  and  widens  (1 0-1 2 cm  wide),  but 
we  see  that  without  a  sloping  RBL  the  plume  is  strongly  deflected 
away  from  the  spreading  axis  such  that  no  plume  material 
reaches  the  ridge  (Fig.  2c). 

Experiment  2  is  performed  with  the  same  plume  source  tem¬ 
perature  but  with  surface  cooling,  which,  combined  with  plate 
spreading,  produces  a  characteristic  upper  thermal/rheological 
boundary  layer  that  slopes  towards  the  spreading  axis.  If  we 
arbitrarily  define  the  RBL  as  the  isosurface  along  which  viscosity 
is  twice  that  of  the  ambient  fluid  (I75^Pas“')  then  the  RBL 
thickness  above  the  plume  source  is  0.25  cm  (that  is,  350  Pa  s"' 
at  15.5  ^C)  after  the  plate-driven  flow  has  been  established.  The 
initial  RBL  slope  in  this  experiment  is  then  T,  assuming  zero 
RBL  thickness  at  the  ridge  axis  14  cm  away. 

The  role  of  the  RBL  on  plume-ridge  interaction  is  apparent 
in  the  paths  of  neutrally  buoyant  tracer  beads  (Fig.  3a).  Most 
beads  are  deflected  by  the  plate  flow,  but  two  ( 1  and  2  in  Fig. 
3a)  migrate  toward  the  ridge,  indicating  successful  plume-ridge 
communication.  The  bead  source  is  positioned  on  the  side  of 
the  plume  heater,  away  from  the  ridge  axis.  Because  of  this,  and 
the  fact  that  only  a  percentage  of  the  plume  is  channelled  to  the 
ridge,  a  large  number  of  beads  (3-7)  track  the  fraction  of  plume 
being  deflected  from  the  ridge.  Long-term  sampling  of  plume 
material  by  the  ridge  is  more  clearly  recorded  by  the  axial  tem¬ 
peratures,  which  increase  steadily  with  time  as  hotter  plume  mat¬ 
erial  reaches  the  ridge  (Fig.  3b).  This  plot  shows  axial 
temperatures  still  increasing  at  95  min,  which  is  45  min  (or 
55  Myr)  beyond  the  arrival  of  bead  1  at  the  ridge.  Spatially,  the 
temperature  anomaly  is  centred  on  RTD  5  with  an  axial  width 
of  roughly  10  cm  (350  km).  This  narrow,  confined  axial  anomaly 
indicates  that  rather  than  spreading  radially  along  the  RBL,  the 
ofif-axis  plume  is  channelled  ridgeward  through  a  narrow  conduit 
as  predicted  from  constructional  volcanism*  and  geochemical 
studies".  The  scaled  anomaly  width  of  roughly  350  km  is  the 
approximate  width  of  geochemical  anomalies  along  the  Mid- 
Atlantic  Ridge  (MAR)"  associated  with  the  Ascension  and  Tris¬ 
tan  hotspots,  both  of  which  are  -^400  km  from  the  MAR,  similar 
to  scaled  laboratory  plume-ridge  offsets  of  500  km. 

Also  consistent  with  the  behaviour  of  the  Tristan  system  is 
the  substantial  cooling  of  the  laboratory  plume  as  it  migrates 
from  the  source  to  the  ridge.  The  laboratory  plume  source  tem¬ 
perature  is  35-40  °C  higher  than  the  ambient  temperature.  The 
equivalent  mantle  plume  temperature  excess  is  ~500  °C  using  a 
mantle  a  value  of  3  x  10“^,  but  at  the  ridge  the  temperature 
anomaly  is  only  ~3  °C,  or  a  mantle  equivalent  of  45  °C,  indicat¬ 
ing  that  substantial  conductive  cooling  of  plume  material  occurs 
between  the  source  and  ridge  axis.  Most  cooling  probably  occurs 
along  the  sub-horizontal  plume  conduit  (path  A-B  in  Fig.  \a) 
where  bead  velocities  drop  to  0.5-1  cm  min"’  and  where  the 
plume  is  in  direct  contact  with  the  cold  upper  boundary  layer. 
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Analogously,  the  Tristan  plume  is  expected  to  have  cooled  by 
as  much  as  200  'C  along  its  migration  path  to  the  MAR^. 

In  experiment  3,  the  surface  is  cooled  even  further  to  yield  an 
initial  RBL  slope  of  —3'  between  the  plume  and  ridge.  Here 
the  resulting  axial  temperature  anomaly  is  broader  (-^20  cm  or 
700  km)  than  with  the  1'  slope  (Fig.  3c).  Assuming  that  the 
axial  anomaly  reflects  directly  the  plume  flux  to  the  ridge,  we 
estimate  that  the  ridgeward  plume  buoyancy  flux  (see  Table  I ) 
is  -^10%  of  the  entire  source  flux,  or  ~3  times  that  estimated  for 
experiment  2,  where  RBL  slope  was  a  third  as  large.  These  results 
suggest  that  the  proportion  of  a  plume  that  feeds  a  nearby  ridge 
depends  directly  on  the  amplitude  of  the  RBL  slope,  which  is 
consistent  with  results  from  two-dimensional  numerical  experi¬ 
ments’®.  Because  the  lithosphere  probably  slopes  by  3-10%  even 
larger  proportions  of  plumes  may  feed  ridges  in  the  mantle. 

Finally,  experiment  4  is  identical  to  experiment  2,  except  with 
a  20%  higher  plume-ambient-temperature  contrast.  Although 
intuitively  it  might  seem  that  increasing  plume  buoyancy  should 
increase  plume  flow  along  the  sloping  RBL  to  the  ridge,  the 
absence  of  any  plume  signal  at  the  ridge  for  this  stronger  plume 
case  (Fig.  3c)  indicates  that  this  relationship  does  not  necessarily 
hold.  Instead,  the  hotter  plume  erodes  a  pocket  in  the  RBL  (Fig. 
la)  which  traps  the  plume  head  in  the  translating  viscous  plate, 
a  result  consistent  with  two-dimensional  numerical  experiments 
(C.K.,  J.-G.  Schilling  and  C.G.,  manuscript  submitted).  The 
trailing  conduit  is  tilted  away  from  the  ridge,  thereby  prohibiting 
further  ridgeward  flow.  Such  behaviour  may,  in  part,  explain 
why  only  the  weakest  plumes  discussed  by  Sleep*’  have  recogniz¬ 
able  signatures  at  nearby  ridges’*.  The  most  robust  present-day 
plume  is  Hawaii,  but  it  is  located  >5,000  km  from  the  nearest 
ridge.  The  prediction  from  these  experimental  results  is  that 
lithospheric  erosion  may  have  prohibited  communication  of  the 
Hawaiian  plume  with  the  East  Pacific  Rise  throughout  its  exist¬ 
ence,  even  during  50-70  Myr  ago  when  it  was  closer  to  the  ridge 
than  it  is  today. 

Our  fully  three-dimensional  variable-viscosity  laboratory  ex¬ 
periments  indicate  that  the  preservation  of  an  upper  rheological 
boundary  that  thins  toward  the  ridge  axis  is  the  primary  require¬ 
ment  for  communication  between  an  off-axis  mantle  plume  and 
a  nearby  spreading  centre.  The  greater  the  RBL  slope,  the  more 
effectively  it  diverts  the  buoyant  low-viscosity  plume  material  to 
the  ridge.  But  extremely  hot  plumes  may  erode  the  RBL,  thus 
precluding  the  plume  from  feeding  the  ridge.  Our  laboratory 
results  show  that  the  plume  signal  at  the  ridge  is  narrow,  indi¬ 
cating  that  communication  could  be  through  a  confined  con¬ 
duit  such  as  has  been  proposed  as  a  result  of  previous 
observations*'".  Continuing  laboratory  experiments  are  investi¬ 
gating  the  effects  of  larger  RBL  slopes  and  wider  ranges  in  plume 
buoyancies  and  plate  velocities  on  plume-ridge  dynamics.  It  is 
important  to  note  that  in  cases  where  ridges  have  migrated  away 
from  ridge-centred  plumes,  the  dynamics  of  interaction  may  be 
affected  by  the  track  of  thin  lithosphere  (e.g.  thermal  groove'*) 
that  is  left  behind.  To  address  this  issue,  the  laboratory  appar¬ 
atus  is  also  being  employed  in  a  study  of  stationary  plumes 
interacting  with  migrating  ridges  and  plumes  beneath  ridge- 
transform  offsets.  □ 
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